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Computer Simulations of High Energy Physics 

Philip John Stephens 

Abstract 

This thesis describes the development of two independent computer programs, Herwig++ 
and Effective. Both of these programs are used for phenomenological predictions of high 
energy physics. The former is used to simulate events as measured at particle colliders. 

The latter is used to generate the mass spectrum of supersymmetric models. 

Simulation of collider events requires the implementation of several different aspects 
of particle phenomenology. After a brief introduction on the relevant aspects of the Stan- 
dard Model and numerical techniques, a new set of variables for parton shower evolution 
are presented. These new variables retain the angular ordering feature of the variables 
found in the original HERWIG software, while improving the Lorentz invariance of the 
shower and improving the coverage of the phase space. These new variables define new 
initial conditions for the shower depending on how the partons are colour connected. 
Also developed is a new model for hadronization. By changing the distribution of prob- 
abilities of cluster decays into hadron pairs this model is able to enforce desired results, 
such as isospin symmetry or meson-baryon ratios, more intuitively. 

The physics of the Herwig-|--|- software is described in detail. The improvements to 
the new evolution variables and the new hadronization model provide are illustrated by 
comparing them against data for e"^e~ events. 

Effective is a program that is able to provide the 1-loop effective potential and 1-loop 
mass matrices for an arbitrary N — 1 supersymmetric model. This software also is able 
to solve the renormalization group equations at one-loop for the parameters of the model 
and, in turn, provide the scale dependent values of these parameters. The program is 
described and some results indicative of its potential are also presented. 
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Chapter 1 



Introduction 

1.1 Field Theory Introduction 

This section briefly introduces a few key ideas tliat are used tlirougliout liigli energy 
pliysics for calculating predictions of physics. I start by introducing the Klein-Gordon 
and Dirac field equations. This is followed by a discussion of the Lagrangian for both 
Abelian and non-Abelian gauge groups. Lastly I explain how calculations are performed 
in perturbation theory. 

We start by introducing the Klein-Gordon field. This field obeys Bose statistics and 
is used to describe all bosons: scalars and vectors. Free fields of this form obey the 
Klein-Gordon equation 

(<9' + m2)0(x) = 0. (1.1) 

The propagator of a field is the Green's function of the field equation. In this case it 
is the Green's function of the Klein-Gordon equation. In order to find this function we 
must introduce a pole prescription for the integral. Figure 11.11 shows the contour used 
for the integral over When the contour is closed below we find the Green's function 
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Figure 1.1: The contour taken to find the retarded Green's function. For xq > yo we 
can close the contour below. For xq < yo we can close the contour above, giving zero. 

only over the range Xq > yo- This is called the retarded Green's function. For x^ > y^ it 



and zero for < 







• W 





Figure 1.2: The contour for the Feynman prescription. When xq > y^ the contour can 
be closed below and when xq < yo the contour can be closed above. 

Using the Feynman prescription of pole contours we find the propagator, known as 
the Feynman propagator, for bosonic fields over all x° and y° is 



IS 




(1.2) 




(1.3) 



This can be written as 




,0 



,0 



(1.4) 
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We now want to see how fields that obey Fermi statistics behave. The Dirac Equation 
is the field equation for spin 1/2 fermions. We first start by introducing the Dirac 
matrices j^. These are in four-dimensional Minkowski space and are given in a popular 
representation by 



^1 0^ 



(T* 



; 7 



-a' 



:i.5) 



where 1 is the 2x2 identity matrix and are the Pauli sigma matrices. Using these 
matrices we define the following notation: ^ = ■j'^a^ and ip = ■^^7°. The Dirac equation 
is then 

{i^ - m)'ijj{x) = 0. (1.6) 

Doing the same as for the Klein-Gordon equation, we can find that the retarded Green's 
function, in momentum space, of a Dirac field is 



Sr^ 



i + m) 

^ — m — rri 



2 ■ 



(1.7) 



Along with the anti-commuting nature of fermions this leads to the Feynman propagator 



Sf{x - y) = 



(27r)^p2 — rn? + ie 



which is 



Sf{x - y) 



-ip- (x-y) 



(0| ^{x)ilj{y) |0) for > y° 



(1.8) 



-{0\i/j{y)ilj{x) |0) for x^ < y° 



:i.9) 



1.1.1 Lagrangian 



In classical mechanics, the fundamental quantity is the action, S. This is the time 
integral of the Lagrangian L. The quantity L can be written as a spatial integration of 
a Lagrangian density, C. It is £ that is generally referred to as the Lagrangian in field 
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theory. This is also done so throughout this thesis. 

The Lagrangian contains the information about all of the fields and interactions of 
a theory. The different types of interactions are mediated by particles which transform 
according to the adjoint representation of a gauge group. In the Standard Model (SM), 
these particles are gauge bosons; spin 1 particles. The kinetic term for a gauge boson, 
A^, is 

-^gauge = ^ ) (-'-■-'-0) 

where 

F^. = - - af'^'^A^^A^^. (1.11) 

In equation ()1.11|) the indices A, B and C are the indices in the adjoint representation of 
the group and the indices /i and v are Lorentz indices, g is known as a coupling constant. 
These define the relative scale of the term in the Lagrangian when compared to other 
terms. The terms f^^'^ are the structure constants of the gauge group. Each group has 
a set of generators t"^. The structure constants for the gauge group are given by the 
relation 

=it^f^^^. (1.12) 

When f^^'^ is equal to for all indices the group is called an Abelian group. When 
they are non-zero it is called a non-Abelian group. From p.ll|) it can be seen that 
for an Abelian group the gauge fields don't interact with themselves. For quantum 
electrodynamics (QED) the group is U{1) which is an Abelian group. The gauge field of 
QED is the photon and therefore we can see that photons don't interact with themselves. 
For quantum chromodynamics (QCD) this group is ^[/(S), which is non-Abelian. It is 
the non-Abelian nature of this group that gives rise to the triplet and quartic gluon 
self-interaction terms which complicate QCD calculations. 

In order for the symmetries of a gauge group to be preserved, all the terms in the 
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Lagrangian must be invariant under the local gauge transformations of that group. The 
Lagrangian is composed of gauge fields, matter fields and derivatives of both. As the 
Lagrangian needs to be invariant under these transformations, only certain combinations 
of these fields are possible. 

The matter fields furnish the fundamental representation of the group. In the Stan- 
dard Model these are fermions or scalars; spin | or spin particles. For a set of local 
transformations, ^^(x), we have the transformation 



Here, the indices i and j indicate the fundamental indices and the index A is the index 
of the adjoint representation of the group. In SU {N) the adjoint representation runs 
from 1 to {N"^ — 1) and N is the number of fundamental indices. In U{N) the adjoint 
indices run from 1 to N'^. For example in SU{2) there are 3 adjoint index values and 
2 fundamental index values. Likewise in SU (3) there are 8 adjoint index values and 3 
fundamental index values. 

We now define the covariant derivative for a local transformation. This is 



Here g is the coupling constant of the group. As mentioned before, these govern the 
relative strength the interactions of one group have when compared with the other 
groups. In order to use the Klein-Gordon and Dirac field equations in the Lagrangian 
we must replace the partial derivate 9^ by the covariant derivative D^. We also require 
the covariant derivative to transform in the same way as the matter fields. This is the 
requirement 




(1.13) 




(1.14) 



A 



(1.15) 
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We can use these transformations to define the transformation of the gauge fields. 
In the Abehan case this is 



A^ix)^A^ix)--d^e{x), (1.16) 
9 

while the non-Abelian case the transformation is more complicated. It is given by 

J^t^A^^nix) (j2t^AAn~\x) + -idM^))n-\x). (i.i7) 

A \ A J ^ 

This requirement prevents the addition of a term like \rn?Afj,A^ to give the field mass. 
Instead this must be done by the Higgs mechanism, which is described in section 11.2.21 

The coupling constants entered into the Lagrangian are known as hare couplings. Due 
to higher order corrections, these are not the same as the physically observed couplings. 
Instead the physical coupling is a renormalized coupling. This means that the higher 
order corrections introduce a shift in the coupling. 

At first it may seem that there are an infinite number of possible terms in the La- 
grangian that can satisfy gauge invariance. It will be shown in Section [1.1.21 that higher 
order terms in perturbation theory will involve integrals over 4-momenta of virtual par- 
ticles. These are divergent integrals. In order to do the calculations these must have 
some cutoff imposed at finite momentum, A. Theories where the observables, expressed 
in terms of the suitably renormalized parameters, have values which are independent on 
A are known as renormalizable theories. Using this it can be shown that theories which 
contain a coupling constant of mass to the negative power are not renormalizable. 

The last constraint on the Lagrangian is that it must have dimension (mass)^. Using 
the constraints of gauge invariance, renormalizability and the correct dimension of the 
Lagrangian the only allowable terms in the Lagrangian, for spinors tp, scalars (p and 
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vectors Afx, are 



All of these terms can have some relevant coupling that has mass dimension larger than 
or equal to 0. For a more thorough discussion of gauge invariance, renormalizability and 
Lagrangian formalism see [1-6]. 

1.1.2 Perturbation Theory 

Earlier we found the Feynman propagator for a field that obeys the Klein-Gordon equa- 
tion. This can be written as 



where T indicates the time-ordering property of the Feynman propagator. This propa- 
gator describes the free field theory. Physical predictions can only be made on theories 
that interact, however. To do so we work in the interaction picture. We now define a 
unitary operator, U{t,to), that takes a field at time to to time t in the presence of an 
interaction. This is known as the interaction picture propagator and is defined as 



where we have divided the Hamiltonian into the free field part, Hq and the interaction 
part ifint- We find U obeys the Schrodinger equation 




(1.18) 




(1.19) 



i^U{t,to) = e^^°(*-*°)(i/ - Ho)e-'"^'-'°^ = Hi{t)U{t,to). 



(1.20) 
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Hj is the interaction Hamiltonian written in the interaction picture. This has the form 



Solving this differential equation with the initial condition U{to, to) = 1 we can find the 
solution as 



where the time-ordering of an exponential is the Taylor series with each term time 
ordered. It is this Taylor series that is used when doing perturbative calculations. Before 
we can define these perturbative calculations we must introduce Wick's Theorem. 

In the interaction picture we can decompose the field into its positive and negative 
energy parts 



(1.21) 




(1.22) 



X 







(1.23) 



The contraction of two fields is defined as 




for x° > 7/°; 




(1.24) 



This is exactly the Feynman propagator that was encountered before 



4>{x)4>{y) = Df{x - y). 



(1.25) 



This allows us to write the time-ordering of fields as 



T {(p{xi)4>{x2) ■ ■ ■ (pi^m)} = N {(p{xi)4>{x2) ■ ■ ■ 4>{xm) + aU possiblc contractions} , 

(1.26) 

where N indicates the normal ordering. This is just the ordering of having all creation 



CHAPTER 1. INTRODUCTION 



9 



operators to the right of all annihilation operators. The identity in ()1.2(ij) is Wick's 
theorem. 

When computing the vacuum expectation value of a time-ordered product any uncon- 
tracted operators from applying Wick's theorem give zero ((0 |A^(any operator)! 0) = 0) 
and the contracted operators are simply Feynman propagators! It is this decomposition 
of the time ordered products that leads to Feynman diagrams. 

When we wish to construct higher order terms we will have states that are created and 
destroyed and never produce observable particles. These are known as virtual particles. 
When computing observables, the contributions due to these particles must be integrated 
over their momenta. As we can see from and ()1.8|1 . these integrals contain for 
each virtual particle and d^p for the integrals over these particles. Using power counting 
it can be shown that some diagrams will contain integrals like 



These integrals need to be bound above*^ by an ultra-violet cutoff, A, in order to be 
computed. Observables must be independent of this cutoff and this condition is what 
leads to renormalizability, as discussed previously. 

The power of perturbation theory is fully exploited when the coupling constants 
are small. This allows us to write ()1.22|) as a Taylor series in order of the coupling 
constant. Using Wick's theorem, we can then decompose the terms of the Taylor series 
into normal ordered products. This allows us to use Feynman diagrams to describe each 
of the normal ordered products. Where perturbation theory is a good approximation 
we only need to evaluate the first few terms of the series in order to approximately 
describe the physics. Unfortunately, the magnitude of the higher order terms cannot be 

'^These also need to be bound below, known as an infrared cutoff. This is another problem that is 
not related to renormalizability. 




(1.27) 
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predicted beforehand. Only by calculating them can one decide how accurate the initial 
calculation really is. 



1.2 The Standard Model 

The Standard Model (SM) is well established as a model that describes the particles 
and all the interactions, except gravity. The predictions of the model have been tested 
to high accuracy by the series of LEP experiments. Recent experimental evidence [7] 
shows that neutrino flavours oscillate which means they must have mass. This is in direct 
contradiction to the SM and is the first evidence of physics beyond the SM. Apart from 
the incorrect description of the neutrino flavour oscillation, the model also predicts the 
existence of the Higgs boson. This has not been experimentally confirmed to date. This 
section explains the particle content of the model, the Higgs mechanism and properties 
of QCD. 

There are theoretical reasons to believe that at higher energies the Standard Model 
will also break down. In section 16.11 1 discuss what shortcomings there are believed to 
be and a theoretical solution to these shortcomings, known as supersymmetry (SUSY). 

1.2.1 Particle Content 

The Standard Model is composed of the SU{3)c, SU{2)l and the U{1)y gauge groups. 
Properties of the interaction eigenstates are governed by these groups. The weak and 
and electromagnetic interactions are not directly governed by SU{2)l and f/(l)y groups, 
however. Instead the mass eigenstates are given by the symmetry breaking of these two 
groups. This will be discussed in section 11.2.21 

The gauge fields for these groups are A^^, WJ^ and B^. The field A is known as 
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the gluon field. As mentioned earlier, these gluons have self-interaction terms. That 
means a gluon, unlike the neutral photon, carries a (colour) charge. These terms make 
calculations with QCD much more complex than with QED. The W and B bosons mix 
through the Higgs mechanism (see section ri.2.2|l to form the W^, bosons and the 
photon, A^. It is the photon that mediates the electromagnetic interaction we observe. 
For each field we have a term in the Lagrangian given by ()1.11|) . 

The matter content of the Standard Model can be summarized quite simply. There 
are leptons, quarks and the Higgs boson. Leptons have SU{2)l and f/(l)y charges and 
the quarks have SU (3)c, SU (2) l and f/(l)y charges. The L subscript of the SU (2) gauge 
means that it only couples to left handed particles. We can see from table iLll that there 
are right and left handed charged leptons, up-type quarks and down-type quarks but 
there are no right handed neutrinos. This means that the left handed leptons interact 
with the W and B bosons, but not with the gluons. A, while the right handed leptons 
only interact with the B boson. The left and right handed quarks interact with the 
gluons and the B boson. The left handed quarks also interact with the W bosons while 
the right handed quarks don't. Table ITTT] shows the charges of the fields in each gauge. 
As will be discussed later, the Higgs boson has SU{2)l and U{1)y charges. 

The strength of the different interactions (QED, Weak, QCD) are dictated by the 
relative size of their coupling constants. In QED this constant is e which is related to 
the SU{2)l coupling, gy/ and the f/(l)y coupling, g' . In weak interactions the SU{2)l 
symmetry is broken and the coupling depends on whether we are coupling to the 
bosons or the boson. Either way this coupling is smaller then e, thus the weak 
interaction is weaker than the electromagnetic one. The coupling constant for QCD is 
g^. This is larger than e; for fields that interact with gluons the QCD terms are most 
often the dominant ones. 

'^Note the different font between the gluon field and the photon field. 
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Field U(l)y SU(2)^ SU(3)^ 

en -1 1 1 

£ -I 2 1 

Q i 2 3 

3 1 3 

dR -I 1 3 

Table 1.1: The fermionic fields of the Standard Model and their charges. There are left 
handed 811(2)^;^ doublets, i and Q, and right handed SU{2)l singlets, e^j, ur and (Ir. 
There is also three families of each type of fermion. 

There is more to this picture than the couplings being constant, however. In fact 
these terms are scale dependent. This is known as the running of the couplings and will 
be discussed in section IT. 2. 31 It is believed that at some large scale, the couplings of all 
the interactions, including gravity, will be of the same size, thus unifying the theories. 

It is known from experiments that these particles have mass. In the SM the fermions 
have right and left handed components which have different charges under the groups. 
Due to gauge invariance, a term like mipLipR cannot be added to the Lagrangian to give 
these fields mass as they have different SU{2)l and f/(l)y charges. Instead, the Higgs 
mechanism is used again to define a Yukawa coupling. This will also be explained in the 
next section. 

1.2.2 Higgs Mechanism 

The Higgs mechanism is a mechanism for generating the masses of the Standard Model 
particles while keeping the Lagrangian gauge invariant. The idea hinges on a scalar field 
being added to the model which has a non-zero vacuum expectation value (VEV). The 
scalar can then couple to the particles in the Standard Model and by doing so defines 
their masses. Because the SM Higgs boson has charges under SU{2)l and U{1)y the 
non-zero VEV breaks the gauge invariance. This is why the B boson of the U (1)y group 
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is not the mediator of the electromagnetic interaction. Instead the photon is, which is 
composed of both the B boson and the W3 boson. 

The SM Higgs boson field, $, is a SU(2)^ doublet and carries U(l)y charge. This 
means that it couples to the W bosons and the B boson through the covariant derivative. 
The field can be defined as 



1.28) 



where U{C,) is a unitary operator in SU{2) with three degrees of freedom. These three 
degrees of freedom are known as the Goldstone bosons. It is these bosons that are 'eaten' 
in order to give the and Z° bosons mass. Though these terms do not appear as 
physical particles they do play a role in calculations depending on the choice of gauge 
fixing term, ^. if is a real scalar field which has a zero VEV and is interpreted as the 
physical Higgs boson field. 

The potential of the Higgs boson is given by a linear combination of <I>^$ term and a 
($t$)2 iQj-m. The result is added to the Lagrangian in the SM. The Lagrangian for the 
Higgs boson is 

>CHiggs = {D^'<^yD^<^ + - A (<l>t$)' . (1.29) 

The positive sign in front of the /x^ is different from the standard mass term in a La- 
grangian. Therefore, when the parameters /i^, A > this gives the potential a minimum 
at ($) 7^ 0. The fi and A parameters also dictate what the VEV of the field is that 
minimizes the potential. The particular combination 



defines the value of the VEV that minimizes the potential. Figure lL3l shows the potential 
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Figure 1.3: The effective potential of tlie Standard Model as v is varied. The three 
lines are the tree level potential, the 1-Loop correction only with the Electroweak par- 
ticles (leptons, W,Z,7,Higgs boson) and the one-loop correction for all Standard Model 
particles. These are all for the combination = 246.0 GeV 

as the VEV is varied for a choice of = 246.0 GeV. It can be seen that at v = 246.0 
GeV this potential is a minimum. Also in the figure is the 1-loop effective potential 
including only the electroweak (EW) particles and the 1-loop effective potential when 
all the SM particles are included. These are also minimized at the same value of the 
VEV. This figure was generated by putting the model into the software Effective, which 
will be explained in Chapter IHl 

The covariant derivative term, Cqov Higgs = (-D^$)*(D^$) from ()1.29|) . when ex- 
panded out is 



Gov Higgs 



( ^ 

\ V-2 J 



:i.3i) 



From this term we find that the boson and the B boson mix. This is known as 
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spontaneous symmetry breaking as the exact symmetries SU{2)j^ and U(l)y are broken 
down to f/(l)em and a broken 5'f/(2)^y. We can see the SU{2)w is broken as the 
bosons and the have different masses. The mixing between W3 and B produces mass 
eigenstates with mass squared | {g'^ + g^r) and 0. These are then interpreted as the 
Z° boson and the photon, A, respectively. This means that the mass eigenstates are 
defined as 



^0^ 
A^ 



^ ^ cos Q^^ — sin Q^^ ^ { Wo^ ^ 



sin Q 



w 



cos 6 



w 



[1.32) 



where cos^ 6w = sin^ 6w 



Here 9w is known as the Weinberg angle. 



The mass squared of the boson is \g'w'^'^) where the bosons is the combination 



V2 



1.33) 



and is the complex conjugate. 

The electric charge, Q, of a field is then given by eQ = eT^ + eY . T3 is the eigenvalue 
of the third generator of the SU{2)l group. For a SU(2)^ doublet this is \ for the 
first component of the doublet and — | for the second component, e is given as e = 
gw^vnOw = g'cosOw We can see from table II. ll that the charged leptons have charge 
—1, neutrinos have charge 0, up-type quarks have charge +1 and down-type quarks have 
charge — |. 

The value of the VEV can be determined from the Fermi constant, Gp- This is 
measured as 1.166 x 10^^ GeV^^. It is defined as 



9w 
8M^ 



1 

2^ 



Gf 
V2 



:i.34) 



and yields the value v = 246 GeV. This value, along with the value of 6w then correctly 
predicts the mass of the boson as well. The fact that this mechanism predicts the 
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masses of the Z'^ and bosons as well as the fact that the photon is massless, in a 
gauge invariant way, has led to the belief that the Higgs boson must exist. Fortunately, 
even though the particle hasn't been found yet its mass is given by niH = Since 
/i is only restricted by loop corrections it gives bounds of the Higgs mass [8] between 
117 GeV and 251 GeV with 95% confidence; experimentally this ceiling has not yet 
been reached. Finding the Higgs boson is one of the main goals of the upcoming Large 
Hadron Colhder (LHC). 

The Higgs boson also allows a gauge invariant way to introduce fermion masses as 
well. Instead of having a term mxipxL'ipxR for each fermionic field x we introduce instead 
a Yukawa matrix Y/j for each fermionic family (leptons and quarks) . The Yukawa term 
in the Lagrangian is 

-^^Yukawa = Y/jtJjjL^^pfR + h.C, (1.35) 

where the subscript L and R denote the left and right handed fields. The VEV of $ 
then generates the mass terms for all the fermions and the mixing between families. Of 
course, we still have independent parameters to define the mass of the fields. 

The most general gauge invariant term that can be added for the quark masses is 

= -Y,jQ^a<PadRj - Y^^e'^'Qia^^Rj + h.C. (1.36) 

The Yukawa matrices are not necessarily symmetric or Hermitian. In fact, there is no 
principle that even requires that they are real valued! However, if CP is conserved, this 
would be true. We can simplify the form of ()1.36|) by diagonalizing the matrices obtained 
from squaring the Yukawa matrices, Y^, Y^. Each one defines two unitary matrices, f/j 
and Wi, for i = u,d hj 

Y,Y^ = U,Dlul Y^Y, = W,Diwl (1.37) 
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where Di is the diagonal matrix with elements which are the positive square roots of the 
eigenvalues. The Yukawa matrices can then be defined as 



[1.38) 



If we now make a chiral rotation of the right handed fields by W and the left handed 
fields by U, these transformations don't affect the couplings of these particles to the 
Higgs field. In this basis, we also find that P, C and T are conserved. 

Since all the up and down type quarks have the identical couplings in QCD the U 
matrices commute with the covariant derivative of QCD. This transformation, however, 
mixes ul and and does affect the SU{2) xU{l) couplings. If we now neglect the 
QCD interactions and write the Lagrangian in the basis of the and ^4, rather 

than the 1^'s and B, we have 



C = e{i^)i + enmen + QmQ + URmuR + dnmdR 



(1.39) 



where 



V2 

1 

7^ 



cos 9 



w 



1 



^L-Y^L + glY ( ~2 ^^^^ dw]eL + CR'y^ sin^ Owen 



+ ul7 



- - sin^ Ow)ul + ur^^" ( sin^ Ow ) ur 

1 1 . O „ \ , T ., / I 



+dLY (^-2^3 ^wjdL + dR-f'' \^ - sin" % ) dR 

2_ 1 - 



:i.4o) 
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These currents contain an abundance of information. For example, the first two show 
that the couple up-type quarks to down-type quarks and neutrinos to charged 
leptons. From the electromagnetic current one can directly read off the electric charges 
of the different particles. 

We can now see how these currents change when the chiral transformations, U and 
W, are applied. We can see that in the electromagnetic current, J^^i, the transformation 
matrices cancel out 

dird^^ ^ (uiy ru^'di = d'^rdi. (1.41) 

This is also true for J^. The current that couples to the W^, however, does change. 
This is 

J^^ = -^<^'di - -^uir {UiUaf d{. (1.42) 
This defines a new matrix 

V = UiUa, (1.43) 

which is known as the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix. This explains 
why strange quarks enter into weak interactions. The boson is able to not only turn 
up-type quarks into down-type quarks but also change the generation in the process. 

The same arguments that were just given can also be applied to the lepton families. 
Since the neutrinos don't interact in any way except by the weak interactions we can 
by convention choose to label the mass eigenstates of the neutrinos according to the 
charged lepton partner it is formed with. Unlike the case of the quarks, there is no 
way to distinguish these states in another way. This way a boson only couples the 
neutrinos of one generation to a charged lepton of the same generation. 
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1.2.3 Running Coupling 

2 

We start by defining a for a theory as a = In QED g is the electric charge of the 
positron and the a is the fine structure constant and is denoted simply by a. In QCD 
g is g^ and the a is labelled as- In order to remove the ultra-violet divergences in the 
perturbative series a renormalization procedure is used. This procedure introduces a 
mass scale /i. Therefore, when we want to calculate a dimensionless physical observable 
at mass scale Q, it can only depend on the ratio Q^/ /i^, which is not constant. The choice 
of 11 is arbitrary and therefore if we were to hold the bare coupling of the Lagrangian 
fixed a physical quantity, R, cannot depend on /x. Instead it must depend only on Q^/n^ 
and the renormalized couplings. We will consider first the running of q;^, as it plays an 
important role in QCD. In QCD this can all be expressed as 



_d_ 
djj/ 



,das d 
dfjA das 



R^O. 



[l.U) 



Identifying 

the /i independence of the observable R for massless particles can be expressed as 



1.45) 



1 + "'"^);^ 



R (e*, as) = 0. 



(1.46) 



It is from this that we define the running coupling as{Q'^) as 



las 



[1A7) 



We can then see that 



dasiQ') 
dt 



dasiQ') OPiasiQ')) 



das 



df3{as) 



:i.48) 
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A solution to ()1.4(ij) is as{Q'^)). Therefore, the scale dependence of R is due solely 
to the running of the coupling. By the same means we can find the running coupling of 
a for QED interactions. 




(a) 



9 

(b) 



Figure 1.4: Graphs which contribute to the QCD (3 function in the one-loop approxima- 
tion (in a physical gauge). 

The running coupling constants are then determined by the renormalization group 
equation (RGE), 

(1.49) 



^2 "^a . 



These always have at least dependence as they are extracted from higher-order loop 
corrections to the bare vertices of a theory. Fig. II. 41 shows contributions to the /3 function 
for QCD at one-loop. Figure IT!^ shows a new interaction due to the non-Abelian nature 
of QCD. The (3 function in QCD then has the perturbative expansion 



p{a) = -bal (1 + b'as + + 0(a|)) . 



:i.50) 



The coefficient b" depends on the renormalization scheme, b and b' don't, however, and 
are given by 



b 
b' 



lie A - 2n 



f 



12n 

11C\ - bCAUf - ^CpUf 
2t:{11Ca - 2nj) ' 



:i.51) 
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where Ca = S, Cf = ^ and ^/ is the number of active flavours. From this we see that 
in order for b to be negative (which corresponds to the leading term of the expansion 
being positive) we must have Uf > ^^y^. For Ca = 3 this is nj > 16. If we look at the 
QED j3 function 

/?QED(a) = T^a' + 0(a3). (1.52) 

OTT 

The leading term of this (3 is always positive. This is where we can see the difference 
the non-Abelian interactions (triplet and quartic gluon vertices) of QCD makes. QCD 
is known as an asymptotically free theory. This means that as becomes smaller as 
increases, corresponding to the opposite sign of the /3 function. 

A fuller discussion of RGBs is given in chapter IHl as a part of the development of the 
software Effective. 



1.2.4 Asymptotic Freedom 

We start first by explaining why the observed charge of the electron decreases with 
distance. Referring back to our previous discussion of the running coupling we see that 
as increases the QED coupling increases. This corresponds to it growing at small 
distances. This can be easily explained as follows. The larger the distance between the 
observer and the charge means that more electron-positron pairs will be temporarily 
created out of the vacuum. These can be considered temporary electric dipoles which 
are preferentially aligned with the positive end towards the electron and the negative 
end away. Fig. 11.51 shows this situation. This effectively screens the bare charge of the 
electron and it appears to have a smaller charge the further away from the electron you 
are. 

In QCD, the picture is not quite so clear. We now have three degrees of colour 
charge, as opposed to one in QED, and non-Abelian interactions. There are two ways 
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Figure 1.5: Virtual e+e pairs are effectively dipoles that screen the bare charge of the 
electron. 

of describing the difference: either as a dielectric effect or as a paramagnetic effect. 

We start by giving the dielectric effect argument. We can define a running charge at 
scale in terms of the charge at an ultra-violet cutoff Auv ^ q^ and a scale-dependent 
dielectric constant, €{q'^). This is 




(1.53) 



and from the previous discussion we find 



1 



= 1 




(1.54) 



This implies that the running charge satisfies the equation 



da{q^) 
dhiq^ 



(1.55) 
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We see that in QED /3 is positive meaning the dielectric constant is greater than one. 
This corresponds to a screening of the charge. In QCD we have the opposite case. Here 
the dielectric constant is smaller than one. This is an antiscreening of the charge. 

If we now assume the vacuum of a relativistic quantum field theory can be treated as 
a polarizable medium we can define a magnetic permeability, /i(g^) that due to Lorentz 
invariance must satisfy the equation 



for all q. We can now analyze the behavour of the magnetic susceptibility, x = ~ 1- 
From ()1.54|) we have 



This can be broken up into two parts. The first term, known as Pauli paramagnetism, 
describes how the spins interact with the magnetic field and the second term, known as 
the Landau diamagnetism, describes how the orbital motion of the particles interacts 
with the magnetic field. When the Pauli paramagnetic term is larger then the Landau 
diamagnetic term the system is considered a paramagnetic system, otherwise it is dia- 
magnetic. In QCD the Pauli paramagnetism term has some dependence on both spin ^ 
quarks and spin 1 gluons. The higher spin gluons make a larger contribution than the 
quarks. The permeability is given by ()1.57|) when 









In QCD a contribution from a particle with spin S" to 6 is given by 




(1.59) 
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This gives a total b term in QCD of 

^_ llCA-2nf 

From this we can see that the permeabihty can only be negative (making it a diamagnetic 
system) when nj > 16. This shows that QCD is asymptotically free due to the colour 
charge carrying spin 1 gluons. 

QCD is an asymptotically free theory. It also has a coupling which runs. The running 
decreases as as increases and therefore high energy QCD allows the methods of 
perturbation theory to be applied. 

This is by no means an exhaustive reference of the Standard Model and more complete 
descriptions are readily available [1-3,9,10]. This section has, however, described the 
basics of the particle content, the Higgs mechanism and two topics of importance in 
QCD: running couplings and asymptotic freedom. The existence of the Higgs boson 
is taken as an assumption for the development of SUSY, and will be done so for the 
remainder of this thesis. 

1.3 Event Generation 

We now move into the discussion of Monte Carlo Event Generators. Chapters El through 
Elare all discussions on the development of the event generator Herwig++. In this section 
I start by introducing the Monte Carlo integration technique and explain why this is so 
useful for high-energy particle physics simulations. I then introduce all the various parts 
of an event: hard process, parton shower, hadronization and hadron decays. 
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The idea of the Monte Carlo approach is that the value of an integral can be calculated 
using random numbers. The same ideas are also applicable to sampling based on a 
distribution. Many calculations of quantum field theory involve matrix elements where 
the amplitude squared is interpreted as a probability; the Monte Carlo approach is 
an excellent fit for computer simulations due to these probability distributions. This is 
because points can be drawn according to a distribution, thus simulating a physical event 
with the correct probabilities. A more thorough discussion of Monte Carlo integration 
is given in [11]. 

1.3.1.1 Simple Monte Carlo Integration 

If we pick random points which are uniformly distributed in a multidimensional 
volume V, the basic theorem of Monte Carlo integration is that the integral of a function 
/ over the multidimensional volume can be approximated by 

jfdVr^V (^{f) ± ■ (1-61) 

The angle brackets indicate the arithmetic mean over the points in V space, 

N N 

in^jpEf'i^^y (1-62) 

i=l i=l 

The "plus-minus" term in (jl.^Hl indicates an estimation of one standard deviation of 
the integral. 

To calculate the value of the integral using p. 611) . one simply generates a set of x 
values in the volume V which are the argument of f{x). Using the two parts of p.62|) 
and the set of values generated the approximate value of the integral is given. 
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There is also an algorithm known as the rejection algorithm which is useful when 
generating values according to a distribution. In this case we don't want to calculate the 
integral but rather generate a point according to the distribution. This can be done by 
simply generating the points uniformly in a volume which is V plus an extra dimension 
bounded by the function. If the point generated lies inside the function we accept the 
point. Otherwise we reject it and generate a new one. If we wanted to take the integral, 
this would then be the percentage of accepted points times the volume sampled over. If 
we consider a function of one variable, /(x), we would generate pairs of points, (x, y). 
We find the ratio of points where y < f{x), and multiply this ratio by the area the 
points are generated in. Figure 11.61 shows an example of this technique. This technique 
is called integration by rejection. In the example above the points are generated in an 
unweighted manner. It is also possible to make the algorithm more efficient or add some 
more information by generating points in a weighted manner as well. 




0.2 0.4 0.6 0.8 1 

Figure 1.6: This shows the randomly distributed points in a square of area 1. The value 
of the integral is the percentage of points under the function 

The method can easily be expanded to integrate over regions with an unknown 
volume. Say we want to integrate over a strange region, C, of which we do not know the 
volume. This can also be easily done with the Monte Carlo technique. Simply expand 
the region the points are generated in to a simpler region, M, which contains C, and 
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the multi-dimensional volume of M is known. When the point lies outside of the desired 
region, set the value of the function to 0. 



1.3.1.2 Non-Uniform Sampling for Monte Carlo Integration 

When working with probability distributions, we will be using the points that lie inside of 
the function to be a physical quantity. Therefore, in the simple approach, there may be 
many points that are generated, that aren't used in the simulation. In particle physics, 
many of this distributions have very sharp peaks and valleys. These peaks and valleys 
cause a large number of points to be generated which aren't used. This is a large waste 
of CPU time and we would like to optimize this. 

Fortunately, this can be done. If we instead generate our points according to some 
other distribution, rather than a uniform one, we can reduce the number of unused 
points. Optimally, we would want to generate them according to the function being 
integrated, but this would require that we know the value of the integral beforehand! 
Instead, we find another function that we can integrate, that is always larger than the 
function we would like to integrate. We want to choose this function so that the sharp 
peaks and valleys of our unknown integral are also, approximately, present in our known 
function. But, in order to sample according to our new function, it must be invertible. 

The evaluation of the integral with this new distribution is not much different than 
for a uniform distribution. Again, we will consider a function of one variable, f{x). Our 
known function is g{x) with integral G. Now we randomly generate the x^s and evaluate 
both f{x) and g{x). A random number, TZ between and 1 is also generated. If TZ is 
is less than the ratio f{x)/g{x), then the point is under the function. Otherwise it is 
not. The value of the integral is then just G times the percentage of the points under 
the function. This technique can greatly reduce the number of unused points that are 
generated. The improvement is dependent on how well the known function, g, estimates 
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f{x). An example of this is shown in Figure ITTfl 




Figure 1.7: This shows the random numbers sampled under a known function, g{x). 
The area of the integral is now G, the integral under g{x), times the percent of points 
under f{x). 

Luckily, there are some general features for QCD matrix elements that can be used 
to take advantage of non- uniform sampling. These are well known [12, 13] and general 
algorithms have been developed to improve the efficiency of generating points for these 
matrix elements. 

There are also some other optimizations that can be performed, such as stratified 
sampling. One program that is particularly useful for particle physics integrals is VEGAS 
[14]. More detail on these optimizations is also given in [11]. 

1.3.2 Hard Process 

There are two separate momentum regimes that different methods of high energy physics 
can be used at. The first is the perturbative regime. This regime is for high momentum 
transfer, short distances. In this regime calculations can be approximated by truncating 
the time-ordered series at any order and applying Wick's theorem. Calculations are 
usually calculated in orders of a (couplings). The hard process fits into this regime. 
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There is also the non-perturbative regime. In this regime calculations are complicated 
to perform because as we saw in Section 11.1.21 calculations depend on a time-ordered 
exponential and each term is itself a complicated expression. As was shown earlier, the 
value of ^5 is largest in the low momentum transfer regime. This means the calculations 
would need to be performed to many orders of as- Currently calculations to next-to- 
next-to leading order (leading order + 2 more orders) are state-of-the-art. In fact, in the 
low momentum transfer regime one of the tools of perturbative physics, the free field 
propagator is not valid anymore. Since non-perturbative calculations are so difficult, 
models of physics are usually used instead of trying exact calculations. The process of 
hadronization is an example of this. 

The hard process is the underlying process that occurs when two beams of particles 
are collided. The description of what happens is divided into two different parts, the 
actual hard process and the Parton Distribution Functions (PDFs). A similar division 
in this discussion is used. Here I discuss the hard process and in the next section explain 
the PDFs and how they affect the hard process. 

The hard process of an event is what describes the interaction of high energy particles 
that 'collide'. There are numerous processes for a given set of incoming particles. For 
example, in an e+e~ collision we could produce another e+e^ pair through a e+e^ 
7* e^e~ or a e'^e^ — Z° — e+e^ process. We could also produce a qq pair in similar 
fashion. In fact, there is an infinite number of possible final states as any number of 
photons or gluons can be emitted by an electrically or colour charged particle. The 
higher order terms fall off rapidly however and contribute only a fraction of the total 
cross section. The actual process that occurs is proportional to the fraction of the total 
cross section for a given process. 

In event generators a specific set of hard processes can be used, without considering 
all possibilities. This way different properties particular to certain hard process can be 
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studied. Given a set of processes, one is chosen based on the ratio of each cross section to 
the total cross sections in the allowable set. At the lowest order of perturbation theory, 
there are relatively few processes to calculate for a given pair of incoming particles. As 
the next order of perturbation will contain a new factor of the coupling, these processes 
are suppressed. It is impossible to know beforehand what the cross section for the next 
order will be; it must be calculated. In particular regions of phase space (soft and 
coUinear) there is a formalism for the splitting of a particle, a, into particles b + c. These 
can then be used to generate approximations to higher order diagrams. More detail on 
this is given in Section 11.3.41 These splitting functions are then used to generate the 
QED and QCD emissions from a parton. This process is known as the parton shower. 

The development of matrix elements for a given process is a complete discussion in 
itself [3,9,10]. Here we only give the results of the main processes studied in this thesis. 
The first is e^e~ ff [4]. The differential cross section for this process with either a 7 
or in the s-channel and 6 the center-of-mass scattering angle of the outgoing fermions 
is 



da 



[(1 + COS0) {Q) - 2QfV,VfFl,{s) 
HAl + V!){A) + Vf)Fi:{s)] 
+ cose ^-AQ J A,AfFlo{s) + %A,V,AjVfFl'{s) 



dcosO 



} 




where 





V2GfMI s{s - M|) 
167ra (s-M|)2 + r|M|' 





s is the centre-of-mass energy and Af and V/ are the vector and axial couplings of 
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Vf = Tj - 2Qfsm^ 9w, Af = Tf. (1.65) 

The functions F].o{s) and -Ffo°(s) are the contributions from the Z^-j interference and 
the Z°-exchange, respectively. This process is of particular importance when comparing 
the simulations to data taken at LEP. For center-of-mass energies well below the 2"° 
peak, the two functions F{s) can be ignored yielding a differential cross section of 



da Tia^Q^ 



d cos 9 2s 
Integrating over 9 gives a total cross section of 



cos^^). (1.66) 



.0 = ^. (1.67) 

Around the Z pole, the F^q{s) term dominates and the cross section is approximately 

z 



To use simulations for hadron-hadron events another set of matrix elements is used. 
Below are the matrix elements squared that have been spin and colour averaged (summed) 
over the initial (final) states for massless partons [4]. 

|2 
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27 sf 



(1.69) 

where ^ is the spin averaged sum and s, t, u are the usual Mandelstam variables for the 
process AB CD: s = {pa +PbY-, i={PA- PcY and -u = (p^ - pof- 

1.3.3 Parton Distribution Functions 

Hadrons are composed of quarks. This means that when two hadrons are collided, it 
is some components of the hadrons that interact fundamentally. Determining which 
part of the hadrons are the ones that interact is a complex issue. This is defined by 
the Parton Distribution Functions. These are developed from the data of deep inelastic 
scattering (DIS) experiments. This section provides a brief discussion of these functions. 
For a more detailed discussion, see [4, 15, 16]. We also discuss how the remaining parts 
of the hadrons are handled in a Monte Carlo event generator. These are called the beam 
remnants. 

The factorization theorem allows the study of the parton constituents to be factor- 
ized into a non-perturbative part and a perturbative part. The non-perturbative part 
is determined from experiments, while the perturbative part can be calculated as a 
perturbation series ordered in the strong coupling constant, as. 

We start by considering the DIS process ep — eX. Figure 11.81 shows this process 
with the relavant momenta k, q and p. We start by introducing the variables 



1^1, .,2 4 (s^ + u" t^ + u^\ 

4^' Igrt 27 iu 3 P ' 
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Figure 1.8: The electron emits a virtual photon which probes the structure of the proton. 

We now define the hadronic tensor, W^". This defines how the hadron will interact with 
the photon. This is given by 

W^'^ip, g) = ^ 5] 5] (p, r 1 7/^(0)1 X) {X | J^(0)| p, r) {2nfS\p + q-Px), (1.71) 

r X 

where is the electromagnetic current and r is the spin of the proton. The variable X 
is summed over all of the possible final products of the process. If we require that this 
tensor conserves parity and we have unpolarized protons, we find that we can decompose 
this unknown tensor into two independent amplitudes, Fi{x,Q^) and F2{x,Q^), called 
structure functions. We now write this tensor as 

W^'^ip, q) = - (^g^'' + F,{x, Q') + I [p' + ^q'^ (p" + ^q"^ F^ix, Q'). 

(1.72) 



CHAPTER 1. INTRODUCTION 



34 



The total differential cross section for this process in terms of these structure functions 
is 



The Bjorken limit is defined as Q^, v oo with x fixed. In this limit the structure 
functions depend, approximately, only on x. This implies the photons scatter off pointlike 
constituents. If we work in the 'infinite momentum frame' we can ignore the mass of the 
proton. Comparing p.73|) to the spin averaged matrix element for the process e~q — > e~q 
yields the relation 



where the hat indicates that the stucture function refers to a quark, not a proton and 
the ^ in this context is the fraction of the protons momentum that the quark constituent 
carries. Measurements show that the momentum is not given by a delta function but 
by a distribution. This means that the quarks carry a range of momentum fractions. 

The ideas generated by studying the structure functions in the Bjorken limit are 
incorporated into the 'naive parton model'. In this model the following assumptions are 
made: 

• fiiO^'^ is the probability that quark of flavour i carries a momentum fraction 
between ^ and ^ + rf^, 

• the photon scatters incoherently off the quark constituents. 
fll.74|l can then be written as 





F2 = xel6{x - C) = 2a;F] 
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Beyond leading order the 'naive' parton model is broken in QCD by logarithms of Q^. 
In the Bjorken limit the transverse momentum is assumed to be small. The higher order 
contributions show that is is not the case though. A quark is able to emit a gluon with 
large transverse momentum with probability 05^^ at large kr- These contributions 
give terms proportional to aslog{Q'^). It is these terms that break scaling. When 
is large enough, these terms compensate for the small value of as- The approximation 
where these terms are summed to all orders is known as the leading-log approximation 
(LLA). It is this approximation that the parton shower is developed at. 

The breaking of the scaling means that the functions fi{x) gain a scale dependence 
and are instead fi{x,Q'^). These can be written at some renormalization scale fi as 



fi{x,n'^) = f^o{x) + 



srfh'{-"(i)'4-.(f)} 



+ ..., (1.76) 



where k is an infrared cutoff and fg is a new function, similar to but describing a 
gluon rather than a quark. Pqq is a splitting function in the case of a quark emitting a 
gluon before scattering off the virtual photon, whereas the Pqg is a splitting function for 
the case of a gluon splitting into a qq pair and the q scattering off the virtual photon. 
fio{x) is an unmeasurable, bare distribution. This distribution absorbs all of the collinear 
singularities at a 'factorization scale' /i. F2{x, Q"^) can then be written in terms of ()1.76p 
by 



«=<?,<? 



Ten. 7 1^^ + 



:i.77) 



This has absorbed all of the finite contribution, C, into the parton distributions. It is 
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possible to factor out an arbitrary finite term from the distributions which leaves behind 
an additional finite contribution. This depends on the 'factorization scheme'. A common 
choice is the MS scheme in which ()1.77j) is expressed as 



i=q,q 



(1.78) 



The functions Cq{x) and Cg{x) are called the coefficient functions and they depend on 
the factorization and renormalization schemes. These are the perturbative part of the 
PDFs. The functions /j(x, yU^) contain all the non-perturbative physics. 

Though the parton distributions are not derivable in the scope of perturbation theory, 
we can see from fll.77|l that the right side cannot be dependent on /i^. Taking the partial 
derivative of both sides with respect to t = /i^, and ignoring the gluon part for simplicity, 
yields the DGLAP equation [17-20] 



This equation is analogous to the (3 functions which describe the variation of as{t) with 
t. Including all the terms this can more generally be written as a {2nq + 1) -dimensional 
matrix equation in the space of the quarks, antiquarks and gluons. This equation is a 
fundamental equation for the parton distribution functions as well as the parton shower. 

When developing an event simulation, the PDFs will define which parton is drawn 
from the hadron for the interaction. The remaining part of the hadron is known as 
the beam remnant. Since the way that hadrons are formed and stay bound lies in the 
non-perturbative regime of QCD and is therefore not well understood, the only way to 
describe the remnants is through phenomenological models. 
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The beam remnant carries a portion of the beam hadron's momentum and often will 
travel straight down the beam pipe, leading to undetectable physics. This isn't always 
the case, however, and these remnants can lead to some detectable physical results and 
as such are worth studying. There are very simple models, such as the UA5 model [21] 
that can be used. This model is simply a parameterization of the data and does not 
scale to higher energies. These simple models lack the adequate ability to describe all 
the effects of the remnants. Instead there is evidence to suggest that these remnant can 
and do interact again within one event. These models are known as multiple interaction 
models [22,23] and are currently still an interesting point of research. The whole process 
of the beam remnants and their interactions is known as the Underlying Event. 

1.3.4 The Parton Shower 

The hard process will generate a set of final-state particles for a given set of incoming 
particles. As the hard process is only accurate up to a given order, it isn't able to 
describe events with high parton multiplicity. As colliders are able to achieve higher 
energies, these high multiplicity events begin to play a larger role in the events. There 
is a need to generate these higher multiplicity final states to some approximation. The 
Parton Shower is able to generate these states in the soft and coUinear regions of phase 
space to all orders. It is also in these soft and coUinear regions which the higher order 
matrix elements are enhanced. 

In this section I explain the parton shower and how it describes the enhanced soft 
and coUinear regions [4]. The parton shower is an approximate perturbative treatment 
of QCD at momentum transfer-squared t greater than some infra-red cutoff tg. This 
treatment is then easy to integrate with a hadronization model which begins at the 
cutoff, and turns the partonic final state into hadronic states. 
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In order to generate the higher multiphcity states we must spht a parton a into two 
partons b + c. We start by assuming 



plpi<^pi = t. (1.80) 

If we have a as an outgoing parton, this corresponds to a timelike shower {t > 0). If 
it is an incoming parton, this is a spacelike shower. For the following introduction we 
will consider the timelike shower. It can be shown that the spacelike shower retains the 
same formulation, it just requires different kinematics. 

We can define the energy fraction as 

z = E,/Ea = 1 - E,/Ea. (1.81) 

Now in the region of small angles, which is the region where the matrix element is 
enhanced, we have 

t = 2EbEc{l - cos 9) = z{l - z)El9^, (1.82) 
and using transverse momentum conservation 

e = = ^. (1.83) 

We then find that the matrix element squared for n + 1 partons, from a g gg splitting 
in the small angle approximation, is 

|Atn+l|' ~ ^CaF{z- e„, 66, €c) \Mn\^ , (1-84) 

where Ca = 3 comes from j^bcjAbc ^j^g function F is the helicity dependent 
splittings. If we average this over the incoming gluon spins and sum it over the outgoing 
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gluon spins we get 



Ca {F) = Pgg{z) = C 



1 — z z 



z 1 — z 



+ z(l-z) 



:i.85) 



This is the unregularized gluon splitting function related to the AltareUi-Parisi kernel [3]. 
If we do the same thing for g qq and q qg vfe find the unregularized splitting 
functions 



Tn [z' + (1 - zY] 



1 + z^ 



:i.86) 
:i.87) 



Here Tr = Tr(t^t^)/8 = 1/2 and Cp = Tr{tH^)/3 = 4/3. Therefore we generally find 



|^„+l|'~^Aa|-Mn|'. 



:i.88) 



We then find the new differential cross section for the n + 1 state is 



dan+l = dcTn^dz^Pbaiz) 

t Zn 



:i.89) 



where Pba{z) is the appropriate splitting function to generate the new state. 

The parton shower doesn't just emit one gluon in the soft or collinear region. It is 
capable of multiple branchings. We start by introducing the Sudakov form factor 



A(t) = exp 



to 



'^-^[dz^Piz) 
t' J 27r ^ ^ 



:i.90) 



Doing this we can write the DGLAP evolution equation ()1.79j) as 



t^fix,t)= / — -p(z)f(x z,t) + I ;/ t— A t , 



;i.91) 
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which is equivalent to 

This has the solution 

fix, t) = A(t)/(x, to) + £ y 1^ f f ^^-^^^ 

From this equation, the Sudakov form factor is interpreted as the probability of evolving 
from to to t without branching. The infra-red singularity of the splitting function is still 
not handled, however, so we must impose an upper limit on z such that z < 1 — e(t). The 
branching for values of z above this are interpreted as unresolvable. These are emissions 
of gluons that are so soft they are undetectable. With the cutoff, the form factor is 
interpreted as the probability of not having any resolvable branchings from to to t. 

This idea of the Sudakov form factors easily integrates into the Monte Carlo method. 
From the interpretation of the form factors we find that the probability of evolving from 
ti to t2 without (resolvable) branching is A(t2)/A(ti). For time-like branchings we can 
find the distibution of t2 by solving 

^ = 7^ (194) 
A(t2) ^ ^ 

where 7^ is a random number in the interval [0, 1]. 

When evolving space-like partons the structure of the hadron must be maintained. 
We can see from ()1.92|) that this requires using the ratio A(tj)//(a;, tj). This means we 
must find the distribution of t2 by solving 

/(x,t2)A(t,)-^' ^'-^^^ 



CHAPTER 1. INTRODUCTION 



41 



In either case if the value of t2 is higher than the hard subprocess, Q^, then we have 
reached the no branching condition for that parton. 

In the timehke case, if a branching does occur, we then need to calculate the mo- 
mentum fraction, z. This is done by solving 



where TZ' is another random number in the interval [0, 1]. To construct the momentum 
of the products of the emission we simply need to generate an azimuthal angle in the 
interval [0, 27r]. If polarization correlations are taken into account, then this angle won't 
be uniform. 

The spacelike case isn't very different. In this case we are evolving backwards and 
trying to maintain the structure of the beam particle which is known from the PDF. 
This requires the modification to the distribution of z of 



Another important aspect of these branchings is angular ordering [4]. Coherent 
parton branching shows that each successive emission must lie within a cone with half 
angle given by the previous emission. This effect is also present in QED radiation. In 
this case it is easier to explain the angular ordering effect. Since the charge partner is 
created with angle, 9, an emission outside of this angle would effectly appear to be an 
emission from a chargeless object. In a sense, it can't be determined which particle this 
emission came from, making it an incoherent emission. Therefore a coherent emission 
from a particle must be emitted within the half angle of the previous emission. A similar 
logic can be apphed to QCD radiation, though it is not as straightforward due to the 
three degrees of colour charge but yields the same result. This means that each emission 




(1.96) 




2tt z 



as P{z) 
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is restricted to have a smaller angle then the previous emission. 

In order to have a coherent parton branching, the shower needs to evolve with a 
variable related to 6, rather then one related to the virtuality. This is given by 

In HERWIG [24] the variable used, which ensures angular ordering, is 

i=E^C>'to- (1-99) 



In Pythia [25] an angular ordered variable is not used. Instead, a virtuality ordered 
shower is used and emissions that violate angular ordering are vetoed. Part of my 
research has been into a new set of evolution variables. This new work is discussed in 
Chapter 121 

The shower does have some more complications, however. It is possible to calculate 
matrix elements to higher order. It is desirable to use these matrix elements as they 
contain all of the terms at a particular order, whereas the shower includes some of the 
terms at next-to-leading log (NLL). As these higher order matrix elements are enhanced 
in the soft and coUinear regions it is desirable to use the shower in these regions and the 
matrix element in the hard regions. In order to do so properly, over- and under- counting 
must be prevented. This process is called the matrix element corrections. 

There are two types of corrections, hard and soft. The soft corrections are when 
the shower approximation is corrected so that it is closer to the hard matrix element. 
This means that a gluon from the shower is prevented from occuring as it has already 
been properly considered in the matrix element. The hard corrections are when the first 
gluon emitted is determined to come from the matrix element, rather than the shower. 
The region where the matrix element is used and the region where the shower is used 
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need to match smoothly in order to correctly describe the physics and these corrections 
ensure the smooth matching of the two. 

The matching of the parton shower and the higher order matrix elements is important 
for generating useful results and ensuring that the simulations describe the theory as 
accurately as it can. There are two different problems with matching higher order 
diagrams. The first is in ensuring that all of phase space is properly covered by the 
next-to-leading order matrix element [26]. There is also a need to ensure the shower 
and matrix element match at next-to-leading log (NLL) without double counting [27]. 
The decays of heay^ quarks and other heay^ coloured particles (such as SUSY particles) 
can also involve the parton shower. The showering of these processes and the matrix 
element corrections to them have also been studied [28]. 

1.3.5 Hadronization 

After the parton shower, we are left with a set of partons that are of the same order 
in virtuality as the cutoff on the parton shower, to- At this stage the interactions 
between the partons become heavily influenced by non-perturbative effects in the low 
momentum-transfer, long distance regime. 

A hadronization model must take the partons from virtuality to down into stable and 
unstable hadrons. We would expect that there would be more hadrons created when to 
is larger, as they have a higher virtuality. Likewise we would expect there to be more 
partons at the initial stage of hadronization if to is smaller. 

Since hadronization models aren't that well understood, it turns out that carefully 
letting the parton shower run to lower scales produces better results. This implies 
that though the parton shower only contains perturbative results, these are often still 
more valid than hadronization models. Quarks and gluons are not observed in collider 
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experiments, hadrons are; therefore in the end we must use a hadronization model to 
study coUider physics. 

The simplest model is to assume that each parton produces hadrons independently 
of the other partons. This is known as the Independent Fragmentation Model. This 
was originally designed by Field and Feynman [29] to approximate scaling of energy 
distributions observed in quark jets in e'^e~ events at moderate energy. This model 
takes a quark and pairs it with an anti-quark in a qq pair drawn from the vacuum. This 
forms a meson and the remaining quark fragments the same way. This continues until 
the leftover energy falls below a threshold. In this model the gluon is split into qq pairs. 
The momentum can be distributed in many ways. All the momentum can be given to 
one of the pair so a gluon behaves just like a quark in this model. The momentum 
could also be given by the g ^ qq AltareUi-Parisi sphtting function. This model has 
several problems. As the final partons are supposed to be on mass-shell it can lead to 
momentum conservation problems. Since there are low energy quarks remaining in this 
model, there is also a problem of colour flow. 

Another model is the String Model. This model assumes that two colour connected 
partons have some colour field between them that grows with seperation. It is usually 
assumed to have a uniform energy per length. This amounts to a linear quark confining 
potential. When this "string" between the two quarks contains too much energy it 
breaks and a qq pair fill each side of the break. This continues until each string is 
considered "stable". At this point each part forms hadrons from the flavours that are 
colour connected. Gluons in this model form kinks in the strings because they carry a 
localized energy and momentum. It is the hadronization of these kinked strings that 
generates results that match experiments better than the independent fragmentation 
model. This is the model that is implemented in Pythia [25] . 

The model that is implemented in HERWIG is the the cluster hadronization model. 
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This model tries to cluster quarks together to form hadrons. The cluster model of 
hadronization relies on the colour preconfinement property of parton branching [30]. 
This property implies that pairs of colour-connected neighbouring partons have a mass 
distribution that falls off rapidly at high masses. It then makes sense to cluster these 
partons into colour- singlet clusters that can decay into observable hadrons. 

The gluons in a cluster model must be split in order to form clusters. This is done non- 
perturbatively as we are in the non-perturbative regime. This means that the gluons 
are just split as a two body decay. This allows the colour connected quarks to form 
colour singlet clusters. It is important to note that in order to do the non-perturbative 
splitting, the gluons have to be given an effective mass. The value of this mass dictates 
the available flavours for the splitting. A common value of 750 MeV is used and this 
allows the gluon to split into u and d flavours. A new cluster model [31] is also used in 
the new event generator SHERPA [32]. 

I have developed a new model based on the old HERWIG cluster model. This new 
model and the improvements it provides are discussed in Chapter IHl 

1.3.6 Decays 

The last part to using the Monte Carlo method in collider simulations is the decays of 
the hadrons. There are thousands of decay modes of hadrons and unstable particles 
from the Particle Data Group [33]. Some of these modes are quite rare, others quite 
common. Since there are so many modes, and most of them are for hadrons, calculating 
a distribution for the decay particles is not easy, or even always possible. 

Though exact calculations aren't always possible, there may be simple distributions 
that are better fits to reality based on certain properties, such as parity violation. In 
Herwig-|--|- a few of these simple matrix elements have been included, in much the same 
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way that they were included in HERWIG. The two commonly used decay matrix ele- 
ments are for three-body decays with the decaying particle of mass mo decaying into 
particles of masses mi, m2 and rris. For free particles that decay weakly the free massive 
{V — Ay matrix elements is used. The decay momenta, Qi, are generated in the three 
body phase space with the weight Vy(mo; mi, m2, ms) given by 

X 4:(ml + mi — S)(S — + rrin) 

W mo; mi, m^, mg = ^"'22 W^' l-^OO 

(mo -I- mj — m2 — m^)'^ 

where S is {q2 + QsY and is given by 

S = n{m2 + m^f + (1 - 7e)(mo - mi)^. (1.101) 



For bound particles that decay weakly this same matrix element is used but it has a 
veto placed on it. This veto is to ensure that the decay products are moving away from 
each other fast enough to no longer be bound. 

Decays aren't always confined to happening after hadronization, however. In the 
SM, for example, the top quark will decay before it ever hadronizes. Also the r will 
decay before it leaves the detector. In fact in SUSY there are many more particles that 
decay before they hadronize. Often these particles are fermions and so spin correlations 
can play an important part in the distributions. I present here a method which can be 
integrated into Monte Carlo simulations to provide the full spin correlations to these 
particles [34]. 

This algorithm is difficult to explain in a completely abstract manner. Instead an 
example will be presented here. Consider a 2 ^ n hard subprocess. Here we label the 
incoming particles ai and 02 and the outgoing particles hi to 6„. The momenta of these 
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P«iK'/«.«^-^KiK2;Ai...A„A<:i«-Al...A^. n ^A,a;, (1-102) 



i=l,n 



where p is the spin density matrix of the incoming particles, Ki are the incoming particles 
helicity, Ai is the matrix element of the 2 ^ n process, Aj is the helicity of the bi and 
D\.x'. is the decay matrix for the bi. Initially all the are just Sx^x'. and the spin density 
matrices are given as Pk,ik[ — ^kik[ for unpolarized incoming particles and 



^1(1 + ^3) ^ 

^ 1(1-^3)^ 



(1.103) 



for longitudinally polarized spin 1/2 incoming particles. Here V3 is the component of 
the polarization parallel to the beam axis. 

Next a bj is chosen at random. The spin density for this particle is given by 

Pa,a;. = ^P«i«;P«24-^«i«2;Ai...a„A<«;4;a;...a; n^AiA^, (1.104) 

where the normalization Np is chosen so the trace of the spin density matrix is one. The 
decay mode of this particle is selected based on the branching ratios. This produces 
particles Ci to with helicity Vi. The momentum of these particles is given by the 
matrix element 

PA,A;-MA,;.,.....-Ml;,.i....;^ n ^kvl- (1-105) 

i=l,m 

Another randomly selected decay product, Ck, is chosen from the decay of bj. The spin 
density for this new decay product is 

P->^K = -^Px.yj^>^r,v,...vmMl..,^_„,^ Y[ Di^^,, (1.106) 



Dp 
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where again the normahzation is chosen so the trace of the spin density matrix is one. 

This process of decaying the products continues all the way up the decay chain until 
a stable particle is reached. Once this occurs the decay matrix is fixed as an identity 
matrix (i.e. Once all of the decay products of a particle have been handled the 

decay matrix for the particle is calculated. Returning to our example, assume that the 
particle bj has had all of its decay products generated. Its decay matrix would be 

i=l,m 

where this too is normalized so the trace is one. 

This whole process continues for all bj until all the decay products have had their spin 
density matrices and there decay matrices calculated. Since all the spin information is 
passed up the chain via the spin density matrices and then passed back down the chain 
via the decay matrices, the whole event has the complete spin correlations built into the 
decay products. For a more thorough discussion and examples of the spin correlation 
effects see [34]. 

There are many packages available that perform selected decays using more advanced 
algorithms. These can almost always be implemented into the Monte Carlo simulations 
at this point. One example of this is EvtGen [35]. This package uses decay amplitudes, 
rather than probabilities, so it can correctly generate the angular correlations in a decay 
chain. Many of the decay modes in this package have been developed from experimental 
data and therefore, match data much better than the simple model built into Herwig++. 
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New formalism for QCD parton 
showers 

2.1 Introduction 

The parton shower approximation has become an important component of a wide range 
of comparisons between theory and experiment in particle physics. Calculations of ob- 
servables that are asymptotically insensitive to soft physics are known as infrared safe 
observables. These can be performed in fixed-order perturbation theory, but the result- 
ing final states consist of a few isolated partons, quite different from the multihadron 
final states observed experimentally. One can attempt to identify isolated partons with 
hadronic jets, but then the energy flows within and between jets are not well represented. 

Currently, the only means of connecting few-parton states with the real world is via 
parton showers, which generate high-multiplicity partonic final states in an approxima- 
tion that retains enhanced coUinear and soft contributions to all orders. Such multi- 
parton states can be interfaced to a hadronization model which does not require large 
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momentum transfers in order to produce a realistic hadronic final state. Hadroniza- 
tion and detector corrections to the fixed-order predictions can then be computed, and 
the results have generally been found to be in satisfactory agreement with the data. 
Infrared-sensitive quantities such as hadron spectra and multiplicities have also been 
described successfully using parton showers. This has strengthened the belief that sim- 
ilar techniques can be used to predict new physics signals and backgrounds in future 
experiments. 

This chapter presents a new shower evolution formalism [36], based on an angular 
variable related to transverse momentum [37-40]. The main aim of these new variables 
is to retain the direct angular ordering of the shower while improving the Lorentz in- 
variance of the evolution and simplifying the coverage of phase space, especially in the 
soft region. The old shower variables used in HERWIG used massless splitting functions 
which created an artificial lower bound for the transverse momentum. This created an 
artificial "dead cone" in the emission from heavier quarks in which no emissions could 
lie. By allowing evolution down to zero transverse momentum and the use of mass- 
dependent splitting functions, the new shower variables permit a better treatment of 
heavy quark fragmentation which eliminates these the sharply-defined coUinear "dead 
cones" . 

In the following section the new shower variables and their associated kinematics and 
dynamics are defined. The appropriate argument of the running coupling, the mass- 
dependent parton branching probability, and the shower evolution cutoff are also given. 
The variables are defined slightly differently for initial- and final-state parton branching, 
and depend on the colour connection of the evolving parton, so in subsequent sections the 
various possible configurations of colour flow between initial and final jets are considered. 

The formalism presented here is implemented in the new Monte Carlo event generator 
Herwig-|-+ [41] which is described in detail in Chapter |3 Results for e"'"e~ annihilation 
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and comparisons with LEP data have been presented in a separate pubhcation [42] 
and are also given in Chapter El The formulae in this chapter could also be used to 
construct a matching scheme for next-to-leading order (NLO) QCD calculations and 
Herwig++ parton showers, similar to that developed for HERWIG showers in [43, 44] 
and implemented in the MC@NLO event generator [45]. 



2.2 New variables for parton branching 

As mentioned in Chapter ^ there are two types of shower evolutions. When a parton 
is space-like {t < 0) it is an initial-state parton and is described here as part of the 
initial-state shower. When a parton is time-like (t > 0) it is a final-state parton and 
is described as part of the final-state shower. Each individual splitting of a parton is 
referred to as a branching. The complete branching history of a given parton is called its 
evolution and the collection of all the evolutions of all the final- (initial-) state partons 
is referred to as the final- (initial-) state shower. 

2.2.1 Final-state quark branching 
2.2.1.1 Kinematics 

Consider parton branching in an outgoing (heavy) quark jet. Define the quark momen- 
tum after the ith gluon emission ^ Qi + ki (see figure EH) in the Sudakov basis as 

qi = aip + PiH + q±i (2.1) 

where p is the jet's "parent parton" momentum (p^ = m^, the on-shell quark mass- 
squared), n is a lightlike "backward" 4- vector (n^ = 0), and q±i is the transverse mo- 
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mentum {q'j_- = — q_|_i , q±i ■ p = q±i ■ n = 0). Then 



A = + - "-"^^ . (2.2) 

2aip ■ n 





Figure 2.1: Final-state parton branching. The blob represents the hard subprocess. 
The momentum fraction and relative transverse momentum are now defined as 



Zi = , P±i = q±i - ZiCi^i_^ . (2.3) 



Then we have 

gli = ^ + T^ + ^^. (2.4) 

Zi 1 Zi Zil^i. Zij 

2.2.1.2 Running coupling 

As mentioned in Chapter [T] the QCD coupling constant, is scale dependent. This 
means that we need to decide what the right scale of a branching is. To find the optimal 
argument of ctg, we consider the branching of a quark of virtuality into an on-shell 
quark and an off-shell gluon of virtuality k"^ [46]. From ()2.4j) . the propagator denominator 
is 



q'-m'=^-^m^ + -^ + ^^ = -^\k' + -[p^^ + {l-zfm^]\ . (2.5) 



z 1 — z z(l — z) 1 — z \ z 
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The dispersion relation for the running couphng is supposed to be 



where Ps(^^) is the discontinuity of as(— fc^). The first term on the right-hand side 
comes from cutting through the on-shell gluon, the second from cutting through the 
gluon self-energy. In our case we have k"^ + [p_|_^ + (1 " zym'^]/z in place of k"^ + /i^. We 
are interested in soft gluon resummation — 1) [40] and so we ignore the factor of 
here. Thus the suggested argument of as is p±^ + (1 — zYm?. In practice a minimum 
virtuality is imposed on light quarks and gluons in the parton shower, and therefore the 
actual argument is slightly more complicated (see below). 

2.2.1.3 Evolution variable 

The evolution variable is not simply since this would ignore angular ordering. To 
have angular ordering, for massless parton branching, the evolution variable should be 
p_i_^/[2(l — zy^ = / z(l — z) [47]. For gluon emission by a massive quark we assume 
this generalizes to (g^ — m?)/ z{l — z). To define a resolvable emission we also need to 
introduce a minimum virtuality Qg for gluons and light quarks. Therefore from ()2.5j) 
the evolution variable is 

= + ^ + (2 7) 

^ z^{l - zf z^ z{l -zY ^ ■ ^ 

where n = max(m, Qg). For the argument of the running coupling we use 



(2.8) 
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The /i term allows for massive quarks to evolve down to P_l < (1 — z)m, i.e. inside the 
dead cone [48,49]. 

Angular ordering of the branching gj — > g-j+i is defined by 

Qi+l < ZiQi . (2.9) 

The factor of Zi enters because the angle at each branching is inversely proportional to 
the momentum fraction of the parent. Similarly for branching on the gluon, /cj A^j+i, 
we require 

< (1 - z,)qi . (2.10) 

2.2.1.4 Branching probability 

For the parton branching probability we use the mass-dependent splitting functions 
given in ref. [39]. These are derived in the quasi- collinear limit, in which p_|_^ and m? 
are treated as small (compared io p ■ n) but p_|_^/m^ is not necessarily small. In this 
limit the q qg splitting function is 



PU^,Pa_') = Cf 



1 + z'^ 2z(l - z)m 



2 



p_i_2 + (1 — zYm? 



(2.11) 



It is the second term of this equation which differs from the splitting functions used in 
HERWIG. Note that at p^ = the factor in square brackets is just 1 — z, i.e. the soft 
singularity at z ^ 1 becomes a zero in the collinear direction. The minimum virtuality 
serves only to define a resolvable emission, and therefore we omit it when defining 
the branching probability in terms of the evolution variable fj2.7|) as 



ds dq^ Cp 2/1 \2~2i^^'^ 



dP{q ^qg) = - — — P^g dz = —as[z (1 - z) q 



27r g2 27r ^ ' ' ' q^ 1 - z 



l + z'- 



2m^ 
z(p 



(2.12) 
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In the case of a final-state gluon splitting into a pair of heavy quarks of mass m, the 
quasi-collinear splitting function derived in [39] is 



R 



1 - 2z(l - z) 



p_i_2 + 



(2.13) 



We note that this splitting function is bounded above by its value Tr = | at the phase 
space boundary p_|_ = 0, and below by Tr/2. By analogy with ()2.7p . in this case the 
evolution variable q is related to the virtuality of the gluon or the relative transverse 
momentum of the splitting by 



zl-z) zH\-z] 



(2.14) 



In terms of the variables g, z^ the g qq branching probability then reads 



dP{9^qq) = ^a,[z\l-zfe]% 



1 - 2z(l - z) 



2m^ 



z{l — z)q^ 



dz 



(2.15) 



In the case of gluon splitting into gluons, the branching probability takes the familiar 
form 



C 



~2i dq^ 



dP{g^gg) = ^as[z\l~zY~q']^^, 



1 - z 



+ 



l-z 



+ z{l-z) 



dz 



(2.16) 



Since we introduce a minimum virtuality for gluons, the relationship between the 
evolution variable and the relative transverse momentum for this splitting is as in ()2.14|) 
but with the heavy quark mass m replaced by Qg. Similarly, for gluon splitting to light 
quarks we use ()2.14|) with = max(m, Qg) in place of m. 
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Consider the initial-state (spacelike) branching of a partonic constituent of an incoming 
hadron that undergoes some hard collisions subprocess such as deep inelastic lepton 
scattering. The momenta are defined as in ()2.H) . with the reference vector p along the 
beam direction. In this case the evolution is performed backwards from the hard sub- 
process to the incoming hadron, as shown in figure 12.21 Thus we now define in place of 





Figure 2.2: Initial-state parton branching. The blob represents the hard subprocess. 



Then 



Zi — J P_Li — ~ ^iCt±i 

a,; 



2 ^2 P±i 



(2.17) 



(2.18) 



We assume a massless variable- flavour- number evolution scheme [50, 51] for con- 
stituent parton branching, setting m = and putting all emitted gluons at the minimum 
virtuality, kf = Q"^. The angular evolution variable now relates only to the angle of the 
emitted gluon and therefore we choose 



-2 _ P-Li^ + ^iQl 



(2.19) 
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with ordering condition simply gj+i < g^. Correspondingly, for the argument of the 
running coupling we now use {1 — zYq"^. 

A different type of initial-state branching occurs in the decay of heavy, quasi-stable 
coloured objects like the top quark. Here the momentum of the incoming heavy object 
is fixed and evolution is performed forwards to the hard decay process. In this case we 
cannot neglect the mass of the parton and (j2.19j) becomes 



+ m 



(2.20) 



while the branching probability ()2.12p is replaced by 



dP{q^qg) = -^a,[{l-zYq']^. 



2ti 



q^ 1 — z 



1 + z' 



2zm? 
q2 



(2.21) 



2.2.4 Allowed regions and termination of branching 

The allowed phase space for each branching is given by requiring a real relative transverse 
momentum, pj^^ > Q. In final-state q —>■ qg branching, we have from ()2.8p 



z\l-zff>il-zff,' + zQl. 



(2.22) 



This yields a rather complicated boundary in the (g, z) plane. However, since 



(i-^)V + ^QXi-^)V, z'Ql 



(2.23) 



we see that the phase space lies inside the region 



m Qg 
— < z < 1 - , 
q q 



(2.24) 
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and approaches these hmits for large values of q. The precise phase space can therefore 
be filled efficiently by generating values of z between these limits and rejecting those 
that violate the inequality ()2.22j) . The resulting threshold for q is slightly larger than 
but of the order of m + Qg. 

In gluon splitting, we obtain the allowed phase space range from ()2.14|1 as 



z- < z < z+, z± = ^yl± ^1- ^ j and q > Afx (2.25) 

where fi = m for splitting into heavy quarks, or /i = max(m, Qg) more generally. There- 
fore, analogously to ()2.24|) . the phase space lies within the range 



^<z<l-^. (2.26) 



Schematically, the parton shower corresponds to selecting a sequence of (g^, Zi) values 
by solving the equations 




(2.27) 



where TZi^2 ^ [0, 1] are uniform pseudorandom numbers. Whenever the algorithm selects 
a value of q below the threshold, branching of that parton is terminated. The minimum 
virtuality Qg thus determines the scale at which soft or coUinear parton emission becomes 
unresolvable. In the absence of such a scale one eventually reaches a region where the 
perturbative expression for the running coupling is divergent. 

One may wish to use a parameterization of as at low scales such that as(0) is finite. 
However, a cutoff Qg is still needed to avoid divergence of the q qg and g ^ gg 
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branching probabilities. Alternatively one could consider parameterizing as such that 
as(0) = 0, e.g. 

as(g') = -f as(Qc) for g < , (2.28) 
where Qc > A. Then the total branching probability below Qc is (for massless quarks) 

Pc{q - qg) = C^^^ C 2z{l + z') dz = , (2.29) 

and no explicit cutoff is required, although of course Qc is essentially playing the same 
role. 

After branching has terminated, the outgoing partons are put on mass-shell (or given 
the virtual mass Qg if lighter) and the relative transverse momenta of the branchings in 
the shower are computed. For final-state gluon splitting we have 

Ip^I = ^z^{l-zyq^-fi^, (2.30) 

or else, if the parent is a quark, 

|p±l = y(l^^?(^V^VP^- (2.31) 

The virtualities of the internal lines of the shower can now be computed backwards 
according to ()2.4|) . Finally, the azimuthal directions of the pj^'s can be chosen [34] and 
the full 4-momenta reconstructed using eqs. ()2.H] and ()2.2j) . 

In initial-state constituent parton branching the evolution is "guided" by the parton 
distribution functions (PDFs) of the incoming parent hadron. Since PDFs are often not 
tabulated below some scale Qs > Qo, one may wish to terminate branching whenever 
q < Qs is selected. In that case the incoming parton is assigned virtuality ~ —Q^ 
and the spacelike virtualities of internal lines are then reconstructed back from to Qq 
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using the transverse momenta deduced from ()2.19|1 inserted in ()2.18j) . 

For initial-state branching in the decay of a heavy, quasi-stable coloured object, 
the branching proceeds in the opposite direction but the reconstruction of momenta is 
similar, using ()2.20p instead of ()2.19|) . 

2.2.5 Treatment of colour flows 

The remaining sections of this chapter present a detailed treatment of the colour flows 
which depend on the choice of the "backward" vector n and on which quantities are 
to be held fixed during jet evolution. Normally n should be taken along the colour- 
connected partner of the radiating parton, and the 4-momentum of the colour-connected 
system should be preserved. The upper limits on the evolution variable q for the colour- 
connected jets should be chosen so as to cover the phase space in the soft limit, with 
the best possible approximation to the correct angular distribution. In setting these 
limits we neglect the minimum virtuality Q"^, which is a good approximation at high 
energies. In the remaining sections we consider separately the four cases that the colour 
connection is between two final-state jets, two initial-state (beam) jets, a beam jet and 
a final-state jet, or a decaying heavy parton and a decay-product jet. 



2.3 Final-final colour connection 

Consider the process a ^ b + c where a is a colour singlet and b and c are colour- 
connected. Examples are e~^e~ —>■ qq and W — > qq'. We need to preserve the 4- 
momentum of a and therefore we work in its rest-frame, 



= Q(1,0,0), = ^Q(l + 6-c,0,A) , p, = ^Q(l-6 + c,0,-A), (2.32) 
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A = A(l, b, c) = y/{l + b-cy - 4b = V(l-& + c)2-4c . 



61 



(2.33) 



For emission of a gluon g from b we write 

Qi = OLiPb + I3in + qu (2.34) 

where qj^^ = k^, qj^^ = — k_L, qj^^ — and we choose 

n=^Q(A,0,-A). (2.35) 

Notice that, if c is massive, the ahgnment of n along pc is exact only in a certain class of 
Lorentz frames. However, if we try to use a massive "backward" vector the kinematics 
become too complicated. 

To preserve Po = ?6 + ?c + % we require 

whereas the mass-shell conditions give 

^ = TTT—T ^ 



A(l + 6-c + A) \a 
where k, = k^^/Q^. Our new variables are 





— ba, 






c 






ba^ 


OLc 


K 






bag^ 


ag 





f^c = \7VTl—r^{—-b'^c] (2.37) 
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where from ()2.7|1 we have 

K= {z'^k-b){l- z)\ (2.39) 



and so \fhjK < z <\. From eqs. ()2.36|) - ()2.39|) we find 



1+6-c+A 



(^1 + 6 - c + z(l - ^)/t + - 6 + c - z(1 - 2)/t]2 - 46 



2 a;!, , , 

1 + — c + A z 



l-z 

Otq Ot}} , 

Z 



with the A's given by (j^TI) . 

2.3.1 Phase space variables 

It is convenient to express the phase space in terms of the Dahtz plot variables 



X, = = (1 + 6 - c)«, + AA . (2.41) 

Q 



Substituting from eqs. ()2.37|) and ()2.40|) . we find 



Xc = 1 — b + c — z(l — z)k, 



Xb = (2-Xc)r+ (z-r)v/a;2 -4c (2.42) 



= (2 — Xc)(l — r) — (2; — r) a/x^ — 4c 



where 

2 V 1 + c — Xc 



r = i ( 1 + ) . (2.43) 



The Jacobian factor is thus simply 



(2.44) 
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dP{q ^ qg) = Cf 



as 



dxh dxr 



27r (1 - 5 + c - Xc) ^/xf^^ 



2b 



1 + z' 



1 — z I — b + c — Xc 



(2.45) 



where 



z = r + 



Xb- {2- Xc)r 
^Jxl - 4c 



r being the function of Xc given in ()2.43j) . 



For emission from parton c we write 



(2.46) 



qi = aiPc + Pin + q±i 



(2.47) 



where now we choose 



n = -Q(A,0,A) 



(2.4^ 



Clearly, the region covered and the branching probability will be as for emission from 
parton b, but with Xb and Xc, b and c interchanged. 



2.3.2 Soft gluon region 

For emission from parton b in the soft region 1 — z = e— i^Owe have 

Xc^l-b + c-ek, Xb^l + b- c - eh' (2.49) 

where 

k' = \ + ^{l-b-c-\) . (2.50) 

Since k is an angular variable, we can express it in terms of the angle 9bg between the 
directions of the emitting parton b and the emitted gluon in the rest frame of a. In the 
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soft region we find 

(l + &-c + A)(l + &-c-Acosgfeg) 

Thus k = b at 9bg = and k oo as 6'{,g — vr. 

For soft emission from parton c, the roles of and Xc, b and c are interchanged. To 
cover the whole angular region in the soft limit, we therefore require k < in jet b and 
k < kc in jet c. We also want the slope of the boundaries to match as they approach 
the soft limit, while still covering the divergence. In this case this gives the condition 



ki kf. 



(2.52) 



and hence 



{kk-b){k^-c) = ^{l-b~c + \f . (2.53) 



In particular, the most symmetric choice is 



«b = ^(1 + 6-c + A) , fte = ^(1-6 + c + A) . (2.54) 



The largest region that can be covered by one jet corresponds to the maximal value of 
k allowed in ()2.4U|) for real i.e. for the maximal b jet 



Kb 



4(1 -2V&-6 + C) . (2.55) 



2.3.3 Example: e+e qqg 

Here we have b = c = p, X = v^^^^^ = "^5 the quark velocity in the Born process 
e~^e~ — i> qq. The phase space and the two jet regions for the symmetrical choice ()2.54|) are 
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shown in figure ESI The region D, corresponding to hard non-colhnear gluon emission, is 
not included in either jet and must be filled using the 0{as) matrix element (see below). 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Xq 



Figure 2.3: Phase space for e~^e qqg for rUg = 5 GeV, = m|, with symmetric 
definition of quark and antiquark jets. 



For the maximal quark jet we get from ()2.55|) 

kg = 4(1 - 2Vp) , (2.56) 



as shown in figure 12.41 together with the complementary antiquark jet region given by 
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Figure 2.4: Phase space for e'^e — > qqg for irig = 5 GeV, = m|, with maximal 
region for the quark jet. 

2.3.3.1 Exact matrix element 



The e~^e — > — qqg differential cross section, where V represents a vector current 
such as a virtual photon, is given to first order in as by [52,53] 



1 (Pay OisCf 



CTy dXq dXq 2tI V 

where 
and 



[X, 



+ 2p)2 + (x,- + 2p)2 + C 



V 



2p 



2p 



(l + 2p)(l-x,)(l-x,-) {l-XqY {l-Xq) 



1) J 



(2.57) 



Cv = -8p(l + 2p) 



c^y = ctq (1 + 2p) V 



(2.58) 



(2.59) 



is the Born cross section for heavy quark production by a vector current, ao being the 
massless quark Born cross section. 



In the case of the axial current contribution e^e A ^ qqg, instead of ()2.57|) we 
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have 



1 (Pa A OisCp 



a A dxq dxq 2ti V 



{x, + 2p)2 + (xg- + 2p)2 + Ca 2p 



2p 



V^{1 - Xq){l - Xq) (l-a;q)2 {I - XqY 



, (2.60) 



where 



CA = 2p[{3 + Xgf-19 + Ap] 



(2.61) 



a A being the Born cross section for heavy quark production by the axial current: 



aA = cToW" 



(2.62) 



2.3.3.2 Soft gluon distribution 



In the soft gluon region 1 — z = e —>■ the branching probability ()2.45|1 becomes 



d'^P ^ as '^'^F ftjA 
dxq dxq 2tt ve^ * 



(2.63) 



In this limit, the exact vector and axial current matrix elements, eqs. fl2.57|) and 
()2.6()|1 respectively, give identical distributions: 



1 d'^av 

cry dXq dXq 



1 d'^aA as 2Cf 



OA dxq dxq 27r ve^ 
l-2p p p 



KK' 



(2.64) 



Since from ()2.50p 



K = V + n 



l-V 
l+V 



> V 



(2.65) 
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the parton shower approximation ()2.6H|1 always overestimates the true result in the soft 
limit, and so correction by the rejection method is straightforward. For small values of 
p we have 




10-3 10-2 10-1 10° 



K 

Figure 2.5: The function /(k) giving the gluon angular distribution in the soft limit, for 
m = 5 GeV, = m\. The exact result ()2.64p . solid curve, and shower approximation 
()2.63|) . dashed, are not distinguishable on this scale. 



2.3.3.3 Dead region contribution 

The integral over the dead region may be expressed as 



— f d'ay^^CpF^ik,) (2.67) 
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where Rq parameterizes the boundary of the quark jet. As shown in figure EEl this is 
actually maximal, but still small, at the symmetric point given by ()2.54j) . 




0.1 1 10 

Kq 



Figure 2.6: The function F^{Kq) giving the contribution of the dead region to the cross 
section, for m = 5 GeV, = m|. Solid: vector current. Dashed: axial current. 



Although the integral in ()2.67j) is finite, the integrand diverges as one approaches the 



soft limit Xn 



1 via the narrow "neck" of the dead region in figure 12.31 or 12.41 This 



could cause problems in generating qqg configurations in the dead region in order to 
apply a matrix element correction [54] . To avoid such problems, one can map the region 
Xq.Xq > I into a region whose width vanishes quadratically illustrated 
in figure ITTI The mapping shown is 



Xq > Xq 



Xq ^ Xq 



1 

4 

1 - 2(1 



7 

4 



X. 



q 5 



x„ 



Xn 



(2.68) 



when Xq > Xq > |. Within the mapped region, the integrand then has an extra weight 
factor of 2(1 — x'^) which regularizes the soft divergence. When Xg > Xq > |, and Xq 
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0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 



Xq 



Figure 2.7: The soft region, with jet boundaries (sohd) and mapped region (dashed), for 
m = 5 GeV, Q'^ = m\. 

are interchanged in both the mapping and the weight. 



2.4 Initial-initial colour connection 



Here wc consider the inverse process + c ^ a where a is a colour singlet of invariant 
mass Q and 6, c are beam jets. The kinematics are simple because we take beam jets to 
be massless: in the cm. frame 



p„ = g(l,O,0), p6 = -g(l,0,l), p,= -Q{l,0,-l). (2.69) 



For emission of a gluon g from h we write 



(2.70) 
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where qj^^ = kj_, qj^„ = — k±, Qj^f, = qj^^ = 0. Notice that in this case the recoil 
transverse momentum is taken by the colour singlet a so we cannot preserve its 4- 
momentum. We choose to preserve its mass and rapidity, so that 

= /3„ = VI + « , (2.71) 

where as before k = k^^/Q^. Now we have 

OLa = OLb-Otg , f3a= f3c- f3g , (2.72) 



and our new variables in this case are 



at, 



(2.73) 



Thus we find 



/3„ - + (1 - ^)'^ , 

-y/l+il-Z^K, 
Z 

l + {l-z)k 

Vi + (1 - zfd ■ 



(2.74) 



2.4.1 Phase space variables 



It is convenient to express the kinematics in terms of the "reduced" Mandelstam invari- 
ants: 

s={qb + qc)VQ\ t = {qb-qgf/Q\ u = {q, - qgf /Q^ . (2.75) 
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1 < s < S/Q^ , l-s<t<0, u = l-s-t (2.76) 

where 5* is the beam-beam cm. energy squared. In terms of the shower variables for 
beam jet b, we have 



s = abPc = -[1 + (1 — -z)/^] 5 t = —C(b/3g = —(1 — z)k , u = —(1 — z)s . (2.77) 



Thus curves of constant k in the (s, i) plane are given by 



k(1 — s) , ^ 

t = — 2.78 

K + S ^ ' 



and the Jacobian factor for conversion of the shower variables to the Mandelstam in- 
variants is 

d{z, K) z 

For the other beam jet c we have t ^ u and thus 

s(l — s) , , 

t = ^ ^ . 2.80 

K + S ^ ' 

We see that in order for the jet regions to touch without overlapping in the soft limit 
s — s> 1, t — > 0, we need k < in jet b and k < kc in jet c, where kc = l/Hb- The most 
symmetrical choice shown in figure 12.81 but we can take or k^. as 

large as we like. 
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Figure 2.8: Beam jets (B,C) and dead region (D) in initial-state branching. 
2.4.2 Example: Drell-Yan process 

Consider radiation from the quark in the Drell-Yan process, qq —>■ gZ^. In the laboratory 
frame we have 

qg = {Pxi,0,Pxi) , qg = {Px2,0,-Px2) (2.81) 

where P = ^VS is the beam momentum. If we generated the initial hard process 
qq — > with momentum fractions Xq, Xq and we want to preserve the mass and rapidity 
of the Z^ we require 

Xi = Xqab , X2 = Xqf3c (2.82) 

where «{, and Pc are given by eqs ()2.74|) . 

The branching probability in the parton shower approximation is 

^ = - , (2.83) 

dzdn 277 K 1 — 2; ' 
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which gives a differential cross section (s = sQ^, etc.) 



1 (Pa D{xi)D{x2) Ois ^ s + u 
(To ds dt D{xq)D{xq) 2^ ^ sHu 



[s^ + {s + uf] 



(2.84) 



where ctq is the Born cross section. The functions D{xi) etc. are parton distribution 
functions in the incoming hadrons; these factors take account of the change of kinematics 
discussed above. 

The exact differential cross section for qq gZ^ to order is 



Since = s + t + u and t < 0, we see that the parton shower approximation ()2.84|) 
overestimates the exact expression, becoming exact in the collinear or soft limit t ^ 0. 
Therefore the gluon distribution in the jet regions can be corrected efficiently by the 
rejection method, and the dead region can be filled using the matrix element, as was 
done in [55]. The benefit of the new variables is that the angular distribution of soft 
gluon emission requires no correction, provided the jet regions touch without overlapping 
in the soft region. As shown above, this will be the case if the upper limits on k satisfy 



2.5 Initial- final colour connection 

Consider the process a + b —>■ c where a is a colour singlet and the beam parton b and 
outgoing parton c are colour-connected. An example is deep inelastic scattering, where 
a is a (charged or neutral) virtual gauge boson. We need to preserve the 4-momentum 




(2.85) 
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of a and therefore we work in the Breit frame: 
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= Q(0,0,-1) , pb = ^Q{l + c,0,l + c) , = ^Q(l + c,0,-l + c) , (2.86) 

where = —Q^, pi = 0, and = cQ^. Notice that the beam parton b is always taken 
to be massless, but the outgoing parton c can be massive (e.g. in W'^d c). 



2.5.1 Initial-state branching 

For emission of a gluon g from the incoming parton b we write 

Qi = otiPb + + g±i (2.87) 
where q_|_g = k_L, q_|_^ = 0, q_|_g = — k_L and we choose 

71 = ^Q(l + c,0,-l-c) . (2.88) 
To preserve Pa = Qc + Qg ~ Qb "we now require 

ab - ac - ag = Pc + Pg - Pb = t^— (2.89) 

-L ~n C 

whereas the mass-shell condition is 

aAQ\l + cf = q^,' + q^ (2.90) 

which gives 

1 + c= — + k(— + —] . (2.91) 



ar \ ar a. 



The new variables for emission from the beam jet are as in ()2.7Hj) . Substituting in 
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(jZnH), we find 



a. 



2z 



(l + c+{l- z)k+ v/[l + c + (1 - z)k]^ - Az{l - z)k) 

1 + c V / 



zab — 



1 + c 



(1 - z)ab , (3b = 



1 c+fl-z)^/? 



1 + c z(l + c)«f, — 1 



(1 - z)R 
[1 + cfab 



(2.92) 



2.5.2 Final-state branching 



Next consider emission from tlie outgoing parton c. In tliis case we write 



(2.93) 



To preserve Pa = Qc + Qg — Qb we require 



ac + ag - ab = Pb - Pc - (3g = 'i- 



(2.94) 



wfiereas tlie mass-sliell condition is now 



aipi iQ'' + ml) = q^.^ + - a\vr? 



(2.95) 



Tlie new variables for emission from an outgoing parton are as in eqs. (I2.38p2.39|) with 
h replaced by c: 



OLr 



etc + ttg ' 



C + 



l-z) 



(2.96) 



Thus in this case we find 



ab = , ttc = z , ttg = 1 — z , 
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Pb = -^[l + c + z{l- z)k] , 

Pc = \^[2c + z{1-z)k] , 
l + c 

^ = ^Y^ya-H- (2.97) 



2.5.3 Phase space variables 



In this process the invariant phase space variables are usually taken to be 



2pa -(lb Pa- Qb 



(2.98) 



In terms of the new variables for emission from the beam parton, we have 



2z(l + c+{l- z)k + ^[1 + c + (1 - z)k]'^ - 4z{l - ^)k)(2^.99) 
Zp = (1 + c)/3e = ^ (l - c - (1 - z)k + V[1 + c + (1 - z)k]^ - Az{l - ^)k)(?.100) 



(1 + c)ab 



with the Jacobian 



d(^Xp, Zp) 



In the soft limit ^ = 1 — e we therefore find for the beam jet 



(2.101) 



1 



Xp ~ 



l + c 



1-e 



ecK 



, ~ 1 - - 



€K 



(2.102) 



and 



d{xp, Zp) e 



d{z,K) (l + c) 



(2.103) 



In terms of the variables for emission from the outgoing parton, 



X 



" (l + c)/36 1 + c + z(1-z)k 



, Zp OLq^ z , 



(2.104) 
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Figure 2.9: Beam jet (B), outgoing jet (C) and dead region (D) in initial-final state 
branching: c = 0.25, ki, = 1.25, kc = 1.5. 



so the Jacobian is simply 



and in the soft limit 



d{Xp, Zp) _ ^ 2 
d{z,k) ~ ^ ' " 



(2.105) 



1 + c 



1 + c 



, Zpr^l-e 



(2.106) 



with the Jacobian again given by ()2.103p . For full coverage of phase space in the soft 
limit we require k < Kf, in jet 6 and k < kc va. ]et c, where 



kbii^c - c) = (1 + c) . 



(2.107) 



Thus the most symmetrical choice is Kf, = 1 + c, = 1 + 2c, as shown in figure 12.91 
On the other hand, any larger or smaller combination satisfying ()2.107p is allowed, as 
illustrated in figure ITTUl for = 10. 
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Figure 2.10: Beam jet (B), outgoing jet (C) and dead region (D) in initial-final state 
branching: c = 0.25, kb = 10, = 0.40625. 

2.5.4 Example: deep inelastic scattering 

Consider deep inelastic scattering on a hadron of momentum P^ by exchange of a virtual 
photon of momentum q^. If the contribution to the Born cross section from scattering 
on a quark of momentum fraction Xb = Q^/2P ■ q is represented by ctq (a function of 
and Q^), then the correction due to single gluon emission is given by 



1 (fa Cpas D{xb/xp) 1 + {xp + Zp - 1] 



aodxpdzp 2-K D{xb) Xp{l — Xp){l — Zp) 
In the soft limit Xp, Zp —>■ 1 we have, from eqs. ()2. 10212. 106|) with c = 0, 

(1 - Xp)(l - Zp) ~ e^k 



f2.108) 



(2.109) 



and so 



1 (fa Cpcis 1 



aodxpdzp n e^K 



(2.110) 
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whereas the parton shower approximation gives 

1 d'^a Cpois 1 

r\j 

(Jo dz dk vr ek 

Since the Jacobian factor ()2.10H) or ()2.105p in this hmit is simply e, the shower approx- 
imation is exact in the soft hmit. 



(2.111) 



2.5.5 Example: qq tt 

We denote the momenta in this process by Pa + Pb ^ Pc + Pd and the 2-^2 invariants 
by 

s = 2pa-pb , i= -2pa -pc , u = -2pa ■ Pd , (2.112) 

so that s + i+u = 0. Colour flows from q to t and anticolour from q to i. Therefore the 
momentum transfer q = Pa — Pc = Pa ~ Pb carried by a colour singlet and we preserve 
this 4-momentum during showering. 

For emission from the incoming light quark or the outgoing top quark, we work in 
the Breit frame for this system, where 

g = Q(0,O,l), pa = ^Q{l + c,0,l + c) , = ^g(l + c, 0, -1 + c) (2.113) 

with = —i— and c = m'^/Q'^. Then the treatment of sects. 12.5.11 and 12.5.21 can 
be applied directly, with the substitution h ^ a since the emitting system is now (a, c) 
rather than (6, c). However, the phase space variables are no longer those of sect. 12.5.31 
since they involve the momenta of the q and t, which in the frame ()2.113|) take the 
general form 

p, = [igV(l + c)2 + 4K,Q^,-ig(l + c)], 
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Pd = [^QV{l + cy + iK,Q±,lQ{l-c)], 



(2.114) 



where K = Q^^/Q^ is related to the 2 — > 2 invariants: 



AK 



:i + c) 



t = -Q^{l + c), u = -s-t, (2.115) 



and so 



1 + 



(1 + c)t 



(2.116) 



For emission from the incoming hght quark we define as in sect. 12.5.11 



Qi = otiPa + I3in + g_Li 



(2.117) 



for i = a,c,g, where q_|_a = 0, q_|_c = — k_L, qj_g = k±, and n is as in ()2.88|) . Then the 
Oj's and /3j's are given by eqs. ()2.92|1 with the substitution 6 — > a. The hght antiquark 
and the antitop are not affected and therefore Qb = Pb, Qd = Pa- This allows the complete 
kinematics of the 2^3 process to be reconstructed. The 2 — >■ 3 invariants can be 
defined as in ref. [44]: 



s = 2qa-qb, ti = -2qa-qc, t2 = -2qb-qd, Ui = -2qa-qd, U2 = -2qb-qc. (2.11^ 



It is convenient to express n = Pc — cq so that (for i = a,c, g) 



q. = («. - cPi)pa + (1 + c)PiPc + q±i 



(2.119) 



Then we find 



S = , ti = aa/?c(l + c)t , t2 = t , Ui = ttaU , U2 = j3c{u — ct) — a^s — 2kx ■ ■ 

(2.120) 
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For emission from the outgoing top we use the results of sect. I2.5.2| again with the 
substitution b —>■ a. Thus we now have for i = a,c,g 



Qi = aiPc + PiPa + q±i (2.121) 
where the a^'s and /Jj's are given by eqs. ()2.97|) with b ^ a, and we find that 

s = (3aS , ti = ac(3ai, t2 = i, ui = (3aU, U2 = OcU - PcS - 2k± ■ . (2.122) 

Similar formulae to eqs. ()2.12()j) and ()2.122|1 . with the replacements a —>■ b and c —>■ d, 
will hold for the case of gluon emission from the colour- connected (qt) system. Using 
these relations, one can study the distribution of gluon radiation in the parton shower 
approximation and compare it with the exact qq tig matrix element. Agreement will 
be good in the soft and/or coUinear regions but there will be regions of hard, wide-angle 
gluon emission in which matrix element corrections should be applied. Alternatively, the 
above equations can be used to formulate a modified subtraction scheme for combining 
fixed-order and parton-shower results, as was done in ref. [44] for a different parton- 
shower algorithm. 



2.6 Decay colour connection 

Consider the process b ca where a is a colour singlet and the decaying parton b and 
outgoing parton c are colour-connected. Examples are bottom quark decay, b cW*, 
and top decay, t bW. Here we have to preserve the 4-momentum of the decaying 
parton b and therefore we work in its rest frame, 

Pfe = mb(l,O,0) , pc = ^mh{l-a + c,0,X) , p<j = ^mb(l + a - c, 0, -A) , (2.123) 
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where a = rri^/ml, c = ml/ml and now 

A = A(l, a, c) = V(l + a - c)2 -4:a= ^/{l-a + c)^ - 4c . (2.124) 
2.6.1 Initial-state branching 

For emission of a gluon g from the decaying parton b we write 

Qi = oiiPb + Pin + g±j (2. 125) 

where q_|_g = k_L, q_|_^ = — k_L, qj_{, = and we choose 

n= ^mb(l,0,l) , (2.126) 
i.e. ahgned along pc in the rest frame of h. The mass-shell conditions give 



/3a = "a, /3c = ttc , f^g = "3, (2.127) 

^5 



with /s: = 'k±^/ml. From momentum conservation 



aa + a^ + ag = — + + — = 1. (2.128) 






Recall that in initial-state branching of a heavy object our new evolution variable is 
given by ()2.2()j) . so we have 

ag = l-z, K = {R-1){1- zf (2.129) 
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where k = /ml > 1. Introducing for brevity the notation 
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w = 1 — {1 — z){k — 1) , M=l + a — c — (1 — z)k , V = \J — Aawz , (2.130) 
from (j2.128|) we find 

u + V 

aa=— , ttc = I - aa - ag = Z - tta . (2.131) 

2w 

2.6.2 Final-state branching 

For radiation from the outgoing parton c we write 

Qi = aiPc + PiU + q±i (2.132) 

where pc is given by ()2.123|1 . Since the colour-connected parton b is at rest in our working 
frame of reference, the choice of the hght-hke vector n in this case is somewhat arbitrary. 
By analogy with the cases treated earlier, we choose it to be opposite to that used for 
the radiation from b, i.e. along the direction of the colour singlet a: 

n = ^mb{X,0,-\) . (2.133) 

The kinematics are then identical with those for final-final connection (sect. I2.3|l . with 
the replacement b ^ c, c ^ a. 
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2.6.3 Phase space variables 



As in sect. I2.3| it is convenient to use the Dalitz plot variables, which in this case are 



Xi 



2gi ■ Pb 



mt 



(2.134) 



For emission from the decaying parton b we have Xi = 2ai + f3i and hence, from ()2.13H) . 



U + V u — V 

Xa = — \- — , Xc = W + Z-Xa, Xg = 2 - w - z = [1 - z)k , (2.135) 

2w 2z 



with the Jacobian factor 



d{z, k) 



(l-z) 



u + V u ~ V ^ a{w — z^ 



2w^ 



2z^ 



vwz 



(2.136) 



In the soft limit 2; — > 1 — e we find 



Xa^l + a — c — ek'f^ , x„ ~ ekh 



(2.137) 



where 



Kb 



A + — (1-a + c-A) . 



(2.138) 



For emission from the outgoing parton c we have, from ()2.42j) with the replacement 
6 — >• c, c — >• a: 



Xa = 1 + a — c — z(l — z)k 



Xc = (2 — Xa)r + (2; — r) a/x^ — 4a 

Xg = (2 - Xa)il -r)-{z-r)^/xl-4:a 



(2.139) 



CHAPTER 2. NEW FORMALISM FOR QCD PARTON SHOWERS 86 
where 

r = -(l + ). (2.140) 

In the soft hmit we have from ()2.50|) 

Xa^l + a — c — eke , Xg ~ ek'e (2.141) 



where 



k;^ = A + ^(l-a + c-A) . (2.142) 



For full coverage of the soft region we require 



which gives in this case 



(2.143) 



(^^5 - l)(/5e - c) = i(l - a + c + Xf . (2.144) 



Note that, while there is no upper limit on Kb, the largest value that can be chosen for 
kc is given by the equivalent of ()2.55j) . 

ke< A{1 + a-2y/c- c) . (2.145) 



2.6.4 Example: top decay 

In the decay t Wbg we have a = m^/m^ = 0.213 and c = ml/mf = 0.026, so for 
simplicity we neglect c. Then for radiation from the top we have from ()2.135|) 



(2.146) 
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where u, v, w are given by eqs. ()2.1H()j) with c ^ 0. The phase space is the region 



0<x„<l — a, 1 — Xn 



1- Xr 



< xw < I + a . 



(2.147) 



Notice that for real xw we require > 4awz, i.e. 



1 < K < 1 + a 





1 


— a) 


a 


1 


-z) 



(2.U8) 



and 



1 — a 



< z <1 . 



(2.149) 



K 



Thus there is no upper hmit on k, but the range of z becomes more hmited as k increases. 
For radiation from the b we have from ()2.139|) 



Xw = 1 + a — z{l — z)k 



Xg = -{2-xw)~{z~-)^Jxl.-4:a. 



(2.150) 



To cover the soft region we require k < Ht for emission from the top quark and k, < Kb 
for that from the bottom, where ()2.144j) gives 



Kb 



- 1 



(2.151) 



The most symmetrical choice would therefore appear to be Kb = Kt — 1 = 1 — a = 0.787, 
as illustrated in figure EHU 

As mentioned above, there is no upper limit on kf. Thus the region covered by gluon 
emission from the top quark can be as large as we like. However, ()2.145j) tells us that 
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Figure 2.11: Phase space for decay t — > Wbg, with symmetric choice of emission regions 
for the b (B) and the t (Ti,T2), and the dead region (D). 

the upper hmit for radiation from the b is 

Rb < 4(1 - V^f = 1.16 , (2.152) 

and correspondingly kf > 1 + ^(1 + y/aY = 1.53. Figure 1^. 1 21 shows this maximal region 
that can be covered by emission from the b, together with the complementary regions of 
emission from the t. 

We note from figs. 12.111 and 12.121 that, for any value of kf, the region for emission 
from the top quark consists of two distinct parts that touch at the point Xg = 1 — ^/a, 
Xw = 2^/0, where the W boson is at rest: a subregion Ti which includes the soft limit 
Xg ^ and a hard gluon region T2. 

The exact t — ^ Wbg differential decay rate to first order in ctg is given in [56]: 

1 d^r _ as Cf f {l + a-xw){l-Xg)+xl 

Fq dxy/dxg IT (1 + a — xw)x'^ | ^ 1 — a 
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Figure 2.12: Phase space for decay t — > Wbg, with maximal region (B) for emission 
from the b, together with complementary regions of emission from the t (Ti,T2) and the 
dead region (D). 



{xw + Xg - 1 - a)"^ 2a(l + a - xw)x 



2(1 -a)2 



+ 



(l-a)2(l + 2a) 



(2.153) 



where Fq is the lowest-order decay rate. In the soft region xw — > 1 + a, — > this 
becomes 

(2.154) 



1 d^T OsCf 

Fq dXwdXg TT Xg 



1 + a — Xw (1 — a)x„ 



For soft gluon emission (1— 2; = e — > 0) from the top quark we have from eqs. (j2.137p2.138j) 



Xw ~ 1 + a — e(l — a) , x„ eh 



(2.155) 



and so the exact form of the soft gluon distribution is 



1 d'T asCp^.^. 
1 dXgdxw TT e 



(2.156) 
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where 

In the same region the parton shower approximation ()2.12|) gives 

d'^P 1 d^P Cp / 1\ asCp^.^. 

1 - - ) = — — fW ■ (2.158) 



dxgdxw {l — a)edzdK it (1 — a)e^K 



Thus we see that, for emission from the top quark, the shower approximation is exact 
in the soft hmit. At higher gluon energies, inside the region Ti the parton shower 
overestimates the exact matrix element and can therefore be corrected easily by the 
rejection method. In the hard gluon region T2, which contributes only a small finite 
correction to the cross section, the parton shower overestimates the matrix element at 
lower values of Xg but underestimates it at the highest values. Therefore a combination 
of rejection and matrix element correction is needed in this region. 

For emission from the bottom quark in the soft limit, we use the results of sect. 12.3.21 
with the substitution b ^ 0, c —>■ a to obtain 



xw ~l+a — eK, Xgr^e[l — a + ) . (2.159) 

-L Ct 



Therefore the exact soft gluon distribution in the b jet should be 

1 d^T as Cf 



To dXgdxw TT 



m) (2.160) 



where 

k 



1 — a)/t 



1 + 



1 - a 



(2.161) 



On the other hand the parton shower approximation in this case gives simply 

d'^P 1 d^P a, Cp 



dXgdxw {l — a)edzdK n {l — a)e'^K 



(2.162) 
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Thus the soft gluon distribution in the b jet region is overestimated by a factor of 

r ~ n 2 



which can be corrected by the rejection method. This factor varies from 1 to 5.2 for 
the symmetric choice of the b jet region kf, = 1 — a depicted in figure 12.111 For the 
maximal b jet shown in figure I2.12| it rises to 8.3. Since the shower approximation is 
exact in the soft hmit for emission from the top, one can reduce the amount of soft 
correction required by decreasing the b jet region and increasing that for top emission, 
in accordance with 1)2.1511) . However, for large values of kt the dead region moves near 
to the coUinear singularity at xw = 1 + a and a large hard matrix element correction 
becomes necessary. 



A new formulation of the parton-shower approximation to QCD matrix elements has 
been presented. This formalism offers a number of advantages over previous ones. Direct 
angular ordering of the shower ensures a good emulation of important QCD coherence 
effects, while the connection between the shower variables and the Sudakov-like repre- 
sentation of momenta ()2.1)1 simplifies the kinematics and their relation to phase space 
invariants. The use of mass-dependent splitting functions with the new variables allows 
an accurate description of soft gluon emission from heavy quarks over a wide angular 
region, including the collinear direction. The separation of showering into contributions 
from pairs of colour-connected hard partons permits a general treatment of coherence 
effects, which should be reliable at least to leading order in the number of colours. Since 
the formulation is slightly different for initial- and final-state showering, formulae for all 
colour-connected combinations of incoming and outgoing partons have been given. 
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This new shower formulation is a key element of the event generator Herwig++ [41]. 
Chapter 0] and El describe how this new formalism is implemented and complete results 
for e^e^ annihilation are presented. 



Chapter 3 



Hadronization 



Hadronization is the process in which the perturbative partons (quarks and gluons) 
from the shower enter the non-perturbative phase and are converted into the observed 
hadrons. This is not well understood and instead is modelled by a hadronization model. 
As discussed earlier there are a few different types of hadronization models. Pythia [25] 
uses a string fragmentation model whereas HERWIG 6.5 [24] uses a cluster hadronization 
model. Here I describe this cluster model and the modifications that have been made 
in the new Herwig++ hadronization model. 



3.1 Cluster Formation 

The shower terminates and leaves colour connected pairs with low virtuahty. These 
partons need to combine to form hadrons. The first step of the cluster hadronization is 
to form clusters out of these colour connected particles. These clusters are made up of 
quark-anti- quark pairs or (anti-)diquark- (anti-)quark pairs. This section presents each 
step of the process of cluster formation up until the decay of the clusters into hadrons. 
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The initial step of the cluster hadronization is to split the gluons into quark-anti-quark 
pairs. Since at the end of the shower the gluons are put on-shell, they must be given 
a mass, rUg, in order to decay into a quark-anti-quark pair. The default value of this 
mass is 0.750 GeV. This is then high enough to isotropically decay into uu and dd pairs. 
Each gluon is randomly assigned a decay into either u quarks or d quarks and is decayed 
uniformly in cos^ and 0. After the isotropic decays the event is left with only colour 
connected (di) quarks and anti-(di) quarks. 

3.1.2 Cluster Formation 

Clusters are themselves colour-singlets. They can be made up of either two partons 
(quark-anti-quark pair) or of three partons (quarks or anti-quarks all of different colours). 
In the cluster model of Herwig-I— |- a three parton cluster will occur in baryon non- 
conserving events, which have colour sources or sinks, and in events where a beam 
particle is a baryon. This can be seen because a cluster composed of three partons has 
baryon number ±1. Since drawing from the vacuum doesn't change the baryon number, 
the decay of these clusters must have baryon number ±1. Unless the cluster is derived 
from a beam baryon, these types of clusters must only occur in baryon non-conserving 
events. 

A cluster is formed simply by finding a colour-singlet pair (or triplet) of partons. 
The 4-momentum of the cluster is just the sum of the momenta of the partons. Clusters 
are created for all sets of colour singlets. If a quark or anti-quark is created from a 
colour source or sink (baryon-violating) it forms a three parton cluster with its colour 
neighbours. Two of these partons are randomly combined into a diquark for the purpose 
of the cluster decay. 
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Primary Light Clusters Primary 6-Clusters 




M/GoV M/GeV 

Figure 3.1: Primary cluster mass distribution in the e^e~ annihilation at various centre- 
of-mass energies, Q, for clusters containing only light quarks (left) and a b quark (right). 

The principle of colour-preconfinement says that the invariant mass distribution of 
the clusters is independent of the centre-of-mass energy. This idea is needed to fully 
separate the perturbative regime from the non-perturbative regime. Figure 13.11 shows 
this distribution for Herwig++. The first plot in figure 13.11 shows the distribution for 
the light clusters. The second plot in figure 13.11 shows the distribution for the b quarks 
only. This shows that the fall off of the distributions is similar for each flavour. 



3.1.3 Cluster Fission 

The mass of a cluster is given by = of the cluster. In order for the shower to 
combine more smoothly with the hadronization, clusters of a large mass are decayed 
into two new clusters. If the shower were to be cut off at a larger scale then there would 
be fewer more energetic partons. On the other hand if it were to be cut off at a smaller 
scale it would produce many lower energy partons. Splitting the clusters with large 
mass into two clusters with smaller mass allows the hadron multiplicity to be much less 
variable with the shower cutoff. In turn this allows the tuning of the shower to not be 
as dependant on the hadron multiplicity, making for more consistent results. 
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Figure 3.2: Cluster Fission. Cluster of mass M decays into two new clusters of masses 
Ml and M2 by drawing a pair from the vacuum. 

A cluster is split into two clusters if the mass does not satisfy the condition 

J\^Clpo„ ^ Cl„,ax'^^^°" + (3.1) 

where Clpo„ and Clmax are parameters and Sc is the sum of the masses of the constituents 
that make up the cluster. 

If a cluster is to be split a qq pair is chosen from the vacuum. Only u,d or s flavours 
are chosen with probabilities given by the parameters Pwtj with i being the flavour. Once 
a pair is chosen from the vacuum the cluster is decayed into two new clusters with one 
of the original partons in each cluster. A schematic of this decay is shown in figure 13.21 
In the case where the partons are not beam remnants, the mass distribution of the new 
clusters are given by 

Ml = mi + (M-mi -m3)ry^, (3.2) 
M2 = m2 + {M - m2 - m^) (3.3) 

where mi and m2 are the two components of the original cluster, is the mass of the 
parton drawn from the vacuum. Mi, M2 are the new cluster masses and are random 
numbers in the interval [0, 1]. P is another parameter of the model. If the parton is not 
a h quark then P is PSpltl otherwise it is PSplt2 (see table 14. Points are only chosen 
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Figure 3.3: The distribution of given in ()3.5|) (blue, dash-dotted) and the actual corre- 
lated distributions (black/red, solid) when ()3.4|1 are applied. This is for a cluster mass 
of 1 GeV, P = 0.5 and all three constituent masses rrii = 0.01 GeV. 

in the range confined by 



M > Mi + Ms; Ml > mi + ms; Ma > ma + mg. 



(3.4) 



The corresponding probability density is 



\P-i 



(M — rrii — TTT's) 



(3.5) 



where N is an arbitrary normalization term. It is important to note that because of the 
constraint in ()3.4|) the distributions of the two clusters are correlated and therefore do 
not exactly follow ()3.5p . An example of this is shown in fig. 13.31 



If parton i is from the beam remnant and there is no underlying event model being 
used, the cluster fissions with a different distribution. This distribution has one parame- 
ter BClpow, which is set to 1.0 GeV by default. This is used to define b = 2/BClpo„ which 
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gives 



rmn 




(3.6) 



where X — M — mi — m2 — 2m^. The mass distribution is then defined by 



/(Mj) = mj + ma 



log^ 
h 



(3.7) 



where 



i?=( 



^rain ~l~ '"l (1 '^f'min)) (j'min ~l~ ^"2 (1 '^f'min)) ■ 



(3.8) 



ri and r2 are two flatly distributed random numbers. This mass distribution is designed 
to decrease rapidly to avoid splitting the cluster many times which would produce large 
transverse energies. This would not be desired as the beam remnant process is meant 
to be a soft process. 



The last stage of the chister hadronization model is the cluster decays. Here the clusters 
of a given flavour (51,52) draw a (di)quark anti-(di) quark pair {q,q) from the vacuum 
and form a pair of hadrons with the combination {qi, q) and {q, 52) ■ The possible hadrons 
are selected based on spin, flavour and phase space. The exact way of accepting and 
rejecting these combinations are described in three variants: HERWIG 6.5 [24] , Kupco's 
method [57] and Herwig++. The method of Herwig++ is described in the next section. 



3.2 Cluster Decays 



3.2.1 HERWIG 6.5 



The method of cluster decays implemented in HERWIG is described here. Before the 
software runs, it initialized a list of data for each type of hadron. This list contains two 
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weights, 

• w: a. spin weight for the hadron. This is where Mgg' is 2s + 1 for the 
largest spin of a flavour group. There are three of these groups which contain aU 
the hadrons of a given flavor combination: (1) uu and dd, (2) ss and (3) all the 
remaining flavours. 

• s: this is a phase space correction weight which is normally set to unity. 

HERWIG 6.5 chooses a mode in several phases. First a flavour is drawn out of the 
vacuum based on the probability 

P, = _ZIH(4_. (3.9) 

The values of Pwt(g) are just parameters for u, d, s, c and b. For the diquarks, however, 
the weight is 

Phi, = ^ (1 + ^^i) Pwt55Pwt,,Pwt,, , (3.10) 

where Pwt^g is the diquark parameter. 

The flavour combinations, (gi, q) and (g, ^2), have the possibility of forming A^^i,^ and 
Nq^q^ different hadrons, respectively. Hadrons a^^^g and bg^q^ are chosen randomly from 
this list. These hadrons have spin weight Wa and Wb and are rejected based on these 
weights. 

Lastly, the resulting pair (a, b) are accepted or rejected based on phase space. The 
probability is given as 

Pphase{a,b) = , (3.11) 

Pmax 

where and are phase space adjusting parameters and p^g_x P* value for 

the lightest hadron pair of the flavours {qi, q) , {q, q2) ■ P* for a two-body decay is the 
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cm. momentum in the decay M — > mim2 and is given by 

, _ y/M^ - (mi + ma) - (mi - ma)^ 

^ " 2M ' ^^-^^^ 

in the region of vahd phase space, M > mi + m2. It is set to zero otherwise. 



In the end all of this gives (approximately) the probability of choosing hadrons a, 



and 6q 52 as 

P{a,„,,b,,,,\q,,q,) = P.J^J^^^ p* " ^^'^^^ 

^''q\,q ^^5,92 ^^^91,9 ^'^^9,92 Pmax 

Technically this isn't the probability because of the rejection scheme used in the algo- 
rithm but it is close and is able to illustrate the main advantages and disadvantages of 
the algorithm. 

The problem with this method, as described by Kupco [57], is that as more hadrons 
are added to the list then the particular flavour content of the new hadrons is suppressed 
by the growing factor A^. In effect, the probability of choosing a hadron of a given flavour 
is proportional to the average of all the p*'s of the flavour. In most majority 
of hadrons are inaccessible due to mass contraints. This leads to a suppression of the 
lighter hadrons of that flavour, even for clusters which are too light to decay into the 
new heavier states. 

To look at this further we analyze the result of isospin symmetric clusters uu, dd, ud 
and du. For a cluster with mass just above threshold for the production of vr^vr^ and 
TT+TT^ these should be produced with the ratio vr*^ : 7r+ : vr^ = 1 : 1 : 1. It is easy to see 
that the ratio between the states generated by ud and du are 1:1. Instead we look here 
at just the part of the ratio which differs from unity, uu and dd. Using ()3.13j) we flnd 
(after assuming Sa = Sb = l and p^o^^o = + = P*max) 

P(7r°,7rV,^) = Pu^^, (3.14) 
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P{n\n-\u,u) = P^^^_^_^, (3.15) 

^^du J-^ud ^^^ud ^^^du 

Pi-n'yM = Pd^_^_, (3.16) 

dd dd 

. p„____. (3.17) 

The M's for uu and rfc? are equal and the M's for ud and rf-u are equal. We can also set 
the flavour probabilities to unity. Lastly, the spin weights of the 7t~^ and vr^ are equal. 
We now find the ratio of pions as 



n ^ 2 W^o 1 1 wi. 



Obviously, these ratios are dominated by the number of hadrons of a particular 
flavour content. It just happens that in HERWIG 6.5 there was the right number of 
MM, dd and ud, du hadrons in the list to give approximately the correct ratios. Figure EiU 
shows the ratio for 7r° : 7r+ as the cluster mass increases for uu and dd clusters. It can 
be seen that the ratios are dictated by the number of uu, dd hadrons that are in the list. 



3.2.2 Kupco Method 

Kupco [57] was the first person to point out the problem of suppression of hadrons when 
new modes of the same flavours were added. He realized that the problem was that the 
probabilities were proportional to an average of the p*'s of a particular flavour content. 

In order to remedy the problem a new set of probabilities for choosing a decay mode 
was used. Instead of splitting the probability into independant parts, as was done in the 
original version, one weight was created for each hadronic mode. 



W{aq^^q,hq^q.^\qi,q2) = PqWaWbSaSbVlfi- 



(3.19) 
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Figure 3.4: The ratio of to 7r+ as the cluster mass is increased and with different 
number of hadrons in the uu list . The vertical bars indicate the threshold for a new set 
of hadrons to be produced. 



The probabihty is then 



W{ag,,g,bq^,q\qi,q2) 



(3.20) 



The addition of new hadrons to the list now increases the probabihty of choosing that 
particular flavour. Because a majority of the hadrons arc quite heavy the p* for any 
mode with the new hadron is zero in many cases. Therefore, this new hadron has no 
effect on the choice of mode for hghter clusters. 

Again, lets look at the example of decaying uu and dd clusters with a cluster mass 
just above threshold for a tt+tt" decay. 



W{n^,TT^\uu) = Puwlo sIopIo , 

W{n',7T'\dd) = Pdwlosloplo^^o, 



(3.21) 
(3.22) 
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W{'K'^,Tr~\uu) = PdW^+w^-s^+s^-pl+^^-, (3.23) 

W{7r+,7r~\dd) = PuW^+w^-s^+s^-pl+^^-. (3.24) 

With the same assumptions as before (sj = l^P^o ^^o — P*+ ^-,Pu — Pd) we get 



P(7r°,7r°|MM) = P(7rO,7rO|rfJ) = , , (3.25) 

P(n+,n-\uu) = P(n+,n-\dd) = _ /g^e) 



2 



Finally this gives us the ratio 



1 1 

7r° : 7r+ : 7r~ = wlo : -w^+w^- : -w^+w^-, (3.27) 



and as the WiS are static a priori weights, we can choose them to be w^o = ;^y^wvillv^ 
which will give the correct ratio of 1 : 1 : 1. 

Now if we allow the cluster mass to increase so that the n^rj mode is accessible, we 
get the same weights for the 7r°7r° and the tt+tt" modes and the new tt^t) mode has a 
weight of 

W{7r°,r]\uu) = PuW^oWrjS^oSrjPlo^^, (3.28) 

and similarly for the dd cluster. The sum of the weights, without some common factors 
and using the same assumptions, is 



E = wlo + 2-w^ow^ ^J" + Wt^+w^-. (3.29) 

This is the same for both the uu and dd clusters. Therefore the probabilities are now 

P(7r°,,r») = ^, (3.30) 
P(Al) = ?1!^4^, (3.31) 
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(3.32) 



To find the ratio we need to make sure we count the particles correctly. Doing so we 
find 



(3.33) 



After forcing the right ratio near the tt+tt threshold by setting the Wj's, we now don't 
have the right ratio after the n^r] threshold. 
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Figure 3.5: The ratio of 7r° to vr"*" as the cluster mass is increased. When new hadrons 
become accessible the ratios differ from one, sometimes drastically. 



Figure 13.51 shows the ratio of 



TT ' 



as the cluster mass is increased. After a 



new production threshold is reached the ratio changes. The change can become quite 
dramatic and the average ratio between thresholds differs from unity. 

Though the ratios vary as the cluster mass changes, this seems overall to be an 
improvement from the original cluster decays. We no longer have the problem that 
the addition of new hadrons to the lists causes decreased probability. We now have 
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increased the probability of choosing a flavour content when new hadrons of that content 
are added. In doing so, however, we have created a new problem. In baryon-conserving 
processes, we only have clusters with a quark-antiquark pair, rather then some diquark- 
quark combination. This means two baryons must be created at a time by drawing a 
diquark-antidiquark pair out of the vacuum. Doing so requires a cluster with a high 
mass, which in turn means a new meson added to the list is likely to be available. 
Consider a uu cluster with high mass. If we have 5 hadrons with ud flavour and 2 
baryons with uud flavour we will have only 4 possible combinations of the uud flavour 
baryons but 25 modes of the ud flavour. If all the modes are accessible this makes the 
baryon probabilities lower, as they are all normalized by the sum of all accessible modes. 
If we then add a 6th hadron to the ud list, we now have 36 combinations. If, again, 
these are all accessible that decreases the probability of creating a baryon even more. 
This problem is addressed and solved in the new method implemented in Herwig++ and 
described in the next section. 



3.3 Herwig++ 

The nonperturbative splitting, cluster formation and cluster fission described in section 
13. II has been implemented in Herwig++. The method for decaying clusters is similar to 
the Kupco method described in section 13.2.21 but has been changed to account for the 
lack of baryon production. Results of the hadronization process in Herwig++ are given 
in this section. Also the changes to the cluster decay model are presented here. 

The concept of colour preconfinement says that the mass distribution of the colour 
singlet systems, after the parton shower, are independent of the centre-of-mass energy 
of the hard process. Figure 13.11 shows the cluster mass distribution generated from 
Herwig++ using the default parameter set. It can be seen that this distribution is 
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Figure 3.6: The cluster mass distribution after cluster fission. This is using the default 
parameter set. It can be seen that the heavy clusters no longer exist and are instead 
folded into the lighter part of the distribution. 

indeed independent of Q, the cm. energy of the hard process. 

The cluster fission stage described in section 13.1.31 truncates this distribution by 
breaking apart clusters that are above the threshold in 1)3.111 . Figure ITT)! shows the new 
distribution after this process. The two bumps around 2 GeV and 5.5 GeV are from 
the charm and bottom clusters. This figure is for the default parameter set [Clpow = 
2.0, Clmax = 3.2 GeV). As can be seen, the heavier part of the distribution is folded into 
the lighter part. 

3.3.1 Herwig++ Cluster Decay Algorithm 

The algorithm used in Herwig++ for cluster decays is very similar to the method pro- 
posed by Kupco. The weights of a particular decay mode is given by ()3.19p . As discussed 
previously, the problem with Kupco's method is that there are not enough baryons pro- 
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I— ^ Q T" CI TTl d'i' c^y 
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wiu. iViULiei 


6 


2.3 GeV 


2.3 GeV 


nig 

^^(IGeV) 


0.750 GeV 


0.750 GeV 


0.630882 GeV 


0.630882 GeV 


CI max 


3.2 GeV 


3.0 GeV 


Clpow 


2.0 GeV 


2.0 GeV 


PSpltl 


1 


1 


PSplt2 


0.33 


0.33 


BlLim 


0.0 


0.0 


ClDirl 


1 


1 


ClDir2 


1 


1 


ClSmrl 


0.40 


0.0 


ClSmr2 


0.0 


0.0 


Pwtd 


1.0 


1.0 


Pwt„ 


1.0 


1.0 


Pwt^ 


0.85 


1.0 


Pwtc 


1.0 


1.0 


PT.T-t- , 


1 n 

-L . U 


1 n 

-L .U 


PWtgg 


0.55 


0.65 


Singlet Weight 


1.0 


1.0 


Decuplet Weight 


0.7 


1.0 



Table 3.1: The parameters for Herwig++. The first group are shower parameters, the 
second are all of the hadronization parameters. The meaning of all the parameters is 
given in chapter 0] 

duced. This is because of the quantity of possible decay modes in the meson sector 
compared to the quantity of modes in the baryon sector. To fix this problem a way to 
separate the two sectors was needed. 

The new algorithm separates the meson sector from the baryon sector. The weights 
from ()3.19|1 are still used, but the sum of weights used for normalizing is only summed 
over modes in either the meson or baryon sector. If a cluster mass is high enough to 
decay into the lightest baryon pair then the parameter Pwtgg is used to decide whether 
to use the baryon sector or meson sector. There is a (1 + Pwtqg)^^ probability of using 
the meson sector and a ^^p^'^'' probability of using the baryon sector. This change not 
only allows for more baryons to be created but also gives more direct control of how 
many baryons are produced through the diquark parameter Pwtgg. 
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In this section the resuhs of the different hadronization schemes of Herwig++ for e^e~ 
events at 91.2 GeV are presented. The values of the parameters in Herwig++ used for 
this study are given in table 13.11 This section presents only a few results related to the 
new hadronization method. Full results of Herwig++ including event shapes and jet 
physics are given in chapter El 

Table shows the results of the new cluster hadronization algorithm in comparison 
to the old algorithm. The column labeled 'Old Model' is the result of using the old 
algorithm with the new parton shower in Herwig++. The column Herwig++ is the result 
of using the new cluster hadronization model with the new parton shower variables. The 
last column, labeled Fortran, is the data generated using the Fortran HERWIG program, 
version 6.5. The data in this table are combined and updated from a variety of sources, 
see ref. [58]. 

Neither of the Herwig++ implementations in the table have been tuned but the 
results in the Fortran column are the result of tuning. 10000 events were used for the 
'Old Model' data and the Herwig++ results are for 100000 events. As the models are 
different a new set of parameters are needed. The parameters for this are also shown in 
table im The results given in the last column, Fortran, are taken from [59]. 

As we can see from the multiplicity results Herwig++ is able to obtain multiplicities 
that are as good as, if not better than, the old HERWIG 6.5 results. The of the data 
sets are given in table 13.31 

We also want to make sure that the momentum distributions of the hadrons match 
those from the data. Xp is defined as 

Xp = (3.34) 
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0.001 


0.001439 


0.0009 


/2(1270) 


D,L,0 


0.168 ± 0.021 


0.113 


0.150273 


0.173 


/^(1525) 


D 


0.02 ± 0.008 


0.003 


0.011739 


0.012 


A,D,0 


0.184 ± 0.018 


0.322* 


0.318519* 


0.283* 


L'*(2010)=^ 


A,D,0 


0.182 ± 0.009 


0.168 


0.180251 


0.151* 




A,D,0 


0.473 ± 0.026 


0.625* 


0.570354* 


0.501 




A O 




n 91 R* 


n 1QA77'^* 
u. 1 t 


n 197 







0.096 ± 0.046 


0.082 


0.066209 


0.043 


J/* 


A,D,L,0 


0.00544 ± 0.00029 


0.006 


0.003605* 


0.002* 




D,0 


0.077 ± 0.016 


0.006* 


0.022621* 


0.001* 


^'(3685) 


D,L,0 


0.00229 ± 0.00041 


0.001* 


0.001775 


0.0008* 



Table 3.2: Multiplicities per event at 91.2 GeV. We show results from Herwig++ with the 
implementation of the old cluster hadronization model (Old Model) and the new model 
(Herwig++), and from HERWIG 6.5 shower and hadronization (Fortran). Experiments 
are ALEPH(A), DELPHI(D), L3(L), OPAL(O), MK2(M) and SLD(S). The * indicates a 
prediction that differs from the measured value by more than three standard deviations. 
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Old Model Herwig++ HEKWIG 6.5 

~All Data from TableljiSI = 543.84/35 = 15.54 = 277.16/35 = 7.92 = 490.52/35 = 14.01 

Table 3.3: results for the different cluster decay methods. 

Xp Distribution for Charged Particles 
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Figure 3.7: Plot of Xp for all flavours. Data from OPAL collaboration. 

where Emax is the centre-of-mass energy. Presented here are the results for Xp for all 
charged particles and ■^^'^'^ = In for all charged particles in uds events. Also the 
momentum distributions of vr^, and are shown. All results are compared to data 
from the OPAL collaboration [60,61]. Again we can see that the results from Herwig++ 
are in good agreement with the data. More detailed results are given in chapter 
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Chapter 4 
HerwigH — h 



The new generation of high energy colhders such as the Large Hadron CoUider (LHC) or 
a future hnear colhder (NLC) require new tools for the simulation of signals and back- 
grounds. The widely used event generators HERWIG [62] and PYTHIA [25] underwent 
tremendous development during the LEP era and have reached the limit of reason- 
able maintenance. Therefore these programs (Pythia?) [63] as well as new projects, like 
SHERPA [32], are being completely (re-) developed in the object-oriented programming 
language C-|--t-. 

Chapters 121 and El introduced two new theoretical improvements to the original HER- 
WIG [24] Monte Carlo event generator. In this chapter I present the implementation 
of these two improvements in the new Herwig++ [41] event generator as well as the 
implementation of the various other parts of the event generator. The Herwig-|--|- event 
generator is built on top of ThePEG [64]. ThePEG is an administrative library which 
defines tools and data structures which are commonly used by Monte Carlo event gen- 
erators. ThePEG and Herwig++ also use CLHEP [65]. This is a package that provides 
general HEP functions. 

This chapter will contain some discussions about object-oriented programming. As 
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this thesis is not intended to be a discussion of object-oriented programming and design, 
detailed discussion of these matters is given only as a reference [66]. For the purposes 
of this chapter a few keywords and very simple conceptual definitions of these terms 
are given here. Again, this is not intended as a complete description, just a guide for 
understanding the text of this chapter. 

Classes are the main components of an object-oriented program. A very simplistic 
view of a class is as a function or algorithm (with accompanying data). In a sense a class 
just serves to implement some functionality or apply some algorithm. A class is much 
more complex and diverse than this, but for the current purposes this should suffice. 

A class can be a particular type of another class. For example, a class could define 
how a particle behaves. It could contain all the data as well as some special functions, 
such as boosts, that are useful when working with a particle. During the shower more 
information is needed for a particle, such as the Sudakov basis quantities a and (3 in 
1)2.111 . Rather then completely redefining a ShowerParticle to behave almost identically 
to a regular Particle, except for the new data, we can instead inherit the old Particle 
class into the new ShowerParticle class. This means that all of the original data and 
functionality of the Particle class is also present in the ShowerParticle class, but we can 
now add the new data and functions for the ShowerParticle. This is useful because in 
the code we can just pass all the Particles and ShowerParticles around together as one 
set of Particles. Then if we want to perform special functions on the ShowerParticles we 
can identify which of the Particles is really a ShowerParticle and apply our operations. 

A function in a class can be defined as a virtual function. Classes that inherit a 
class with a virtual function can redefine that function. These functions can actually 
be defined to be non-existent. If this is done the class is known as an abstract class. 
This allows for the definition of an interface without defining the implementation. For 
example, there are many things that all matrix elements must have; all matrix elements 
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must have a set of incoming and outgoing particles. But the actual values of the matrix 
elements, for a given set of incoming and outgoing particles with momenta, are different. 
This is a perfect candidate for abstraction. A class ME can be defined that provides 
the definition of a function double me2 (incoming, outgoing) but does not define it. 
Then any class that inherits ME can implement the function. If we now think about an 
algorithm that uses a matrix element, such as computing a cross section, we can define 
this algorithm independent of what matrix element we want to use. It just needs to 
know that the matrix element class will return a double when given a set of incoming 
and outgoing particles (or momenta). 

This brief discussion of a few terms from object-oriented programming should be 
adequate to comprehend the rest of this thesis. The rest of this chapter describes the 
implementation of the various parts of the Monte Carlo event generator: hard subpro- 
cess, PDFs, parton shower, underlying event, hadronization, and decays. 



4.1 Hard Subprocess 

As discussed in the first chapter, there is only a small set of matrix elements currently 
implemented into Herwig-I— 1-. For e'^e~ annihilation there is just one. It is the e'^e" — > 
7*/Z° — > qq matrix element. Though this is not itself a QCD matrix element it inherits 
the ME2to2QCD class from ThePEG. This is because the ME2to2QCD class has a function 
which returns the number of accessible flavours for a given scale, which is needed for the 
matrix element at a given scale. 

There is a much larger set of matrix elements for pp collisions. All of these inherit 
the ME2to2QCD class. Provided with ThePEG is the set of matrix elements: MEqq2qq, 
MEQG2QG, MEGG2GG, MEQQ2QQ, MEGG2QQ, MEQQ2qq, MEQQ2GG and MEQq2Qq. The lower and 
upper case q are used to decipher processes that have different flavours. Also imple- 
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merited with Herwig++ is a Drell-Yan matrix element qq '-/* / ^ i^i^ . This type 
of process is particularly useful for studying the initial-state shower. 

Implementation of the hard subprocess requires both a matrix element and a phase 
space sampler. The method of sampling phase space is important for the efficiency of the 
Monte Carlo. As discussed before, matrix elements'^ can have several peaks and valleys. 
Using a uniform sampling for this is extremely inefficient and advanced samplers can 
improve the efficiency drastically. The default sampler used in Herwig++ is known as 
the ACDC sampler. This is a component of ThePEG. ACDC is an acronym of Auto- 
Compensating Divide-and-Conquer Phase Space Generator [67]. This algorithm uses a 
divide-and-conquer scheme to divide the phase space into uniform sections which have 
different maxima. Figure 14.11 shows an example of a function which has been divided 
into two sections, each of which is sampled uniformly. 
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Figure 4.1: This figure shows a function being integrated which has been divided into 
two regions, Uq and Ui. 

The way the algorithm works is to generate points uniformly in phase-space. When 
a point is generated in which the matrix element is larger than the value being sampled 

''By 'matrix element' we actually mean 'square modulus of the matrix element' 
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under, ujq, the space is divided into two regions. A region is found in which the function 
is larger than ujq and this region is then sampled under uJi. The auto-compensating part 
of the the algorithm compensates for the fact that the peak has been undersampled and 
oversamples the new region until it is consistent. The shaded region of fig. 14.11 shows 
this new region which is auto-compensated. 



4.2 PDF 

Chapter ^ introduced the concept of the parton distribution functions. As mentioned 
earlier there are different ways to parameterize the non-perturbative component of the 
PDF. Two such implementations are the Gliick-Reya-Vogt (GRV) [68] PDFs and the 
Martin-Roberts-Stirling-Thorne (MRST) [69] PDFs. Both of these have been imple- 
mented in Herwig++. 

The implementation of the PDFs requires the implementation of four virtual func- 
tions inherited from the PDFBase class. These functions are: 

• bool canHandleParticle(tcPDPtr particle) const; this function is used to 
specify what hadrons this PDF can work with. 

• cPDVector partons (tcPDPtr particle) const; this function is used to specify 
what partons can be extracted from the given hadron. 

• ApproxMap approx(tcPDPtr particle, const PDFCuts &) const; this function 
is used to specify, approximately, the upper limits of the parton densities, of each 
parton, for the given hadron. The PDFCuts class is used to give as input the 
kinematical region that the PDF will be used in. 

• double xfKtcPDPtr particle, tcPDPtr parton, Energy2 partonScale, double 
1, Energy2 particleScale = 0.0*GeV2) const; this function is the actual im- 
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plementation of the PDF. This function returns the momentum fraction as I — 
log(l/x). 

Optionally, the PDF can also be specified with the functions: 

• double xfxCtcPDPtr particle, tcPDPtr parton, Energy2 partonScale, 

double X, double eps = 0.0, Energy2 particleScale = 0.0*GeV2) const; 
this function is used to specify the PDF and return the momentum fraction as x. 
By default, this just returns exp (— xf l(. . . )). 

• double xfvKtcPDPtr particle, tcPDPtr parton, Energy2 partonScale, 
double 1, Energy2 particleScale = 0.0*GeV2) const; 

this function returns just the valence part of the PDF with momentum fraction 
I — log (l/x). By default this just returns 0. 

• double xfvx(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale, 
double X, double eps = 0.0, Energy2 particleScale = 0.0*GeV2) const; 

this function is the same as xf vl except it returns the momentum fraction as x 
instead of L It also returns by default. 

The GRV PDFs have been implemented as part of ThePEG. The GRV94L and 
GRV94M PDFs have been implemented. These are two implementations of the GRV 
method where optimal fits of the distributions of the partons have been made to different 
data. 

The MRST PDFs have also been implemented as part of Herwig++. This implemen- 
tation is a wrapper to a previous C++ implementation [70]. The original implementation 
was able to read in the data from a file. The files are standardized and when new data 
is analyzed and integrated into the distribution, new data files are created. The new 
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versions are easily integrated into the event generator by simply reading in these data 
files. 

In order for the event generator to function correctly something must be done with 
the remnant of the incoming hadron left behind when a parton interacts. These are 
handled by RemnantHandlers. ThePEG has a default remnant handler, BaryonRemnant, 
but this is designed to be integrated with the Pythia? string fragmentation model. This 
remnant handler works in two ways. If a valence parton is chosen by the PDF then the 
only remnant is the colour connected diquark remaining in the beam particle. If a sea 
parton is used in the hard subprocess then a colour connected parton is produced along 
with a colourless hadron. The implementation of BaryonRemnant in ThePEG requires a 
Px-generator and a 2;-generator. Though the 2;-generator can be designed to work with 
other hadronization models, it is really a feature of the string fragmentation model and 
isn't something that is needed in Herwig++. 

Instead, since Herwig++ uses a cluster hadronization model, we use a different rem- 
nant handler. This remnant handler, also called BaryonRemnant, is really nothing more 
than a place-holder in Herwig-I— 1-; the initial-state shower in Herwig-I— |- is designed to 
evolve back to a valence parton. This class simply generates a diquark, with either spin 
1 or spin 0, with a flavour depending on the parton drawn by the PDF. This remnant 
is often wrong because the initial-state shower evolves backwards to a valence parton 
but the generation of this valence parton is not known at the time the remnant is cre- 
ated. Correctly handling the remnant is instead implemented as part of the backwards 
evolution and the original remnant created by BaryonRemnant is replaced^ 

''BaryonRemnant is only implemented as it is a requirement of ThePEG. 
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In chapter 01 have developed the theoretical background of the parton shower and how 
it takes the matrix elements calculated perturbatively and evolves down to a stage where 
hadronization models are valid. Chapter |21 presented the development of a set of evolu- 
tion variables used in the Herwig++ shower. In this section I present the implementation 
specific features of the Herwig4-+ shower. This includes how the shower is initialized, 
how it terminates and how the momenta of the particles are set. 

4.3.1 Hard Matrix Element Corrections 

Given an n-jet process it is sometimes possible to calculate the matrix element for n + 1 
jets and this matrix element can be used in the Monte Carlo. There are already soft 
and coUinear jets of order n + 1 from the shower of the n-jet matrix element and these 
jets must be matched with the n + 1-jet matrix element to avoid double counting. This 
matching is known as hard matrix element corrections. 

Hard matrix element corrections have only been implemented for the process e^e~ 
'-f*/Z^ qq. To make a hard correction a pair of three-body phase space variables, x, x 
are generated according to the original two-jet matrix element. Emissions in the dead 
region of fig. 12. HI are only accepted according to the three-jet matrix element. In the 
case of e~^e~ —>■ qq only 3% of the events are corrected by the hard matrix element. 

If a hard emission is added the gg-final state is replaced by a qqg-finaA state and 
the orientation of either the quark or the anti-quark is kept with weights and x^, 
respectively. This results in properly oriented three-jet events, apart from finite mass 
effects [71]. This procedure takes into account the most important subleading higher- 
order corrections that are not given by the parton shower. 
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The ShowerVariables class has a variable, MECorrMode, which sets whether the hard 
matrix elements are on or off. 

4.3.2 Initial Conditions 

As we saw in chapter ^ there is an angular ordering principle that restrict the showering 
to occur only in the cone of half angle between two colour connected partons. Since 
the shower is restricted to a cone in relation to one of its colour partners, the first step 
of the shower is to determine which colour partners the soft and coUinear showers will 
occur between. There is a flag. Approach, used in the PartnerFinder class that sets the 
way the shower partner is chosen. If this flag is zero, then the partner is set completely 
randomly amongst the partons that are colour-connected to it and the partners of all 
partons are set independently of each other. This means if we have particle a, which 
chooses its shower partner to be particle b, particle b does not have to choose a as its 
shower partner. Instead, if it had colour partners a and c it would randomly choose 
between these. On the other hand, if Approach is set to 1 then it randomly selects a 
shower partner and sets both particles to be shower partners with each other. 

A partner is chosen for each gauge 'charge' of the parton. For example, a quark has 
a QED charge and a QCD charge. A colour partner and an electric charge partner are 
both set. This allows the QED showers to compete directly with QCD ones. 

As shown in chapter El the evolution of a particle is carried out in the Sudakov basis. 



where p is the momentum of the particle which is evolving and n is a lightlike (n^ = 0) 
vector with 3-momentum in the 'backwards' direction, which is conventionally set to 
that of the colour partner of the particle in their cm. frame. q± is in the transverse 



q = ap + l3n + q±, 
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direction and satisfies q± ■ p = q± ■ n = 0. 

Once the partner is chosen the initial value of the evolution variable, q, is set. The 
value of this variable depends on whether both partners are initial-state, final-state or 
a combination. If one parton is initial-state and one is final-state the values are 



qi = \/{Pi+Pf) +m}, (4.2) 
qf = ^{p^+p^f + 2m}, (4.3) 

where pi is the momentum of the initial-state parton and Pf is the momentum of the 
final-state parton. This corresponds to ()2.107|) for the most symmetrical choice, where 
q^ = kQ^. For final-final shower partners the initial conditions are 



qi = ViPi+P2)-iPi+n2), (4.4) 
q2 = V{pi+P2) ■ {p2 + ni), (4.5) 

where pi and p2 are the momenta of the incoming partons. rii has three momentum 
equal to pi in the cm. frame and nf = 0; similarly for n2. This corresponds to ()2.54j) 
for the symmetric choice. Lastly, for the choice of two initial-state particles the initial 
conditions are 



qi = q2 = \/iPi+P2f- (4.6) 



This corresponds to the symmetric case of ki = K2 = 1, shown in fig. 



4.3.3 Initial-State Shower 

The initial-state shower evolution begins with the two incoming partons that have been 
chosen from the PDF. These partons are considered as on-shell partons in the hard 
matrix element calculation and the initial q are set as described in the previous section. 
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Once the initial value is set for the evolution variable, each parton then evolves inde- 
pendently of each other. The evolution of one parton proceeds using the veto algorithm. 
For each possible type of branching, a — >■ 6c, a new q and a z are generated based on 
ratios of the Sudakov form factor 



Rba {Qi^Qi+l] 



fbix/z,q)Aba {qcQi 



fh{x,q)Aba {qc,qi+i)' 
with ff, the PDF and Aba the Sudakov form factor 



A,, (qc, q) = exp ^ - £ ^ I dz'^^^^P^a q) 6 (px > 0) [> . 



(4.7) 



(4i 



qc is the lower cutoff and by default is set to be the non-perturbative gluon mass. 



rrin 



750 MeV. The running coupling, as(z, g), depends on the evolution scale and the 



momentum fraction. The argument is (1 — z)'^(f' for reasons given in chapter El The 
Pba are the quasi-coUinear splitting functions for the massive partons [72]. For QCD 
branchings these are 
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(4.9) 
(4.10) 
(4.11) 
(4.12) 



and m^ = for initial-state radiation. For the QED case, we change ag to the fine 
structure constant aem and use the branching for q ^ q'y. Ignoring the parton mass this 



is 



Pliz,q) = el 



1 + z' 



(4.13) 



where e„ is the electric charge of the parton, in units of elementary charge. 
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The Q{p± > 0) function in ()4.8j) is used to ensure that it is possible to reconstruct 
transverse momentum, p±, from the evolution variables, q determines the relative trans- 
verse momentum. For quark branching this is 



Eqn. ()4.7|l gives the probability of no branching above the scale gj+i. 1 — Rl^{qi, qi+i) 
is therefore the probability for the next branching to happen above gi+i. The derivative 
with respect to gj+i is then the probability density for the next branching to happen at 
the scale qi+i. 

Since fl4.7|l is not directly solvable, the veto algorithm is used. In this algorithm each 
part of the distribution is sampled independently by a function which is always greater 
than the desired one, and a veto is placed on the emission if the ratio of the actual 
function to the approximated function is larger then some random number. This gives 
several veto points 




(4.14) 



For the gluon branching, in the initial state shower, this is 




(4.15) 



Wi = e(p_L > 0), W2 




(4.16) 



where 




(4.17) 






(4.19) 
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ggg{z) = Ca + 



1 



(4.20) 



z 



9, 





(4.21) 



1 - 



z 



The trick to the veto algorithm, to ensure that things are correctly sampled, is that 
first a g is generated according to 



with Qs the current scale, and the upper and lower bounds on z, respectively and 
I{zq) the value of Qba integrated over z and from zero to 2:0. A z is then generated 
according to the approximate functions, Qba between z+ and z_. Each veto is then 
tested. If a veto fails a new q is generated. This time, however, instead of in ()4.22j) 
it is equal to the q that was vetoed. If this q is smaller than qc then there is no more 
branchings. The largest value of q generated from each of the branchings decides which 
of the branchings to use. 

Once a splitting has been chosen a new initial state parton and a new final state 
parton are created. The final state parton is taken as an on-shell parton. This will 
be put off-shell during the final-state evolution. The momentum fraction, z, from the 
Sudakov form factor and the new scale q are passed to the new initial state parton so it 
can spht. This process is repeated until no new scale is chosen below qc for any of the 
branchings. From the interpretation of the Sudakov form factor, this means that there 
is no more branching and the evolution of this parton has terminated. 

As discussed before, beam remnants aren't handled properly by the remnant handler. 
Instead they are handled here. When the initial-state evolution of a parton terminates 
a set of forced branchings are imposed. There are three types of termination points: a 
valence quark, a gluon, or a sea quark. If the evolution terminates on a valence quark. 




(4.22) 
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then a diquark^ of the appropriate flavours is produced as the beam remnant. If instead 
the evolution terminates on a gluon, a forced splitting of the gluon from a quark of a 
valence flavour is imposed. The momentum fraction z is distributed according to the 
splitting function. The q is just distributed by dP = dq/q between qc and qs. The 
remnant is again just a (di)quark of the flavour(s) that remain in the beam particle. 
The most complex case is when the evolution terminates on a sea quark. In this case 
two forced splittings are imposed. The first is to force a splitting into a gluon. Again 
the momentum fraction is distributed according to the appropriate splitting function. 
The q is also distributed as dP = dq/q between qc and (jg. The new gluon is then forced 
to split into a valence quark in the same manner described above. 

This process is repeated for both incoming partons. Once they have evolved and 
their remnants are correctly set the momenta of all the partons needs to be set. The 
initial condition is that the beam particles are coming in with known momenta. The 
Sudakov variable a of the first parton is then set to unity, the (3 is set to zero and its 
X is known from the backwards evolution. Each child of the parton then has its a set 
to zai, for the initial state particle and (1 — z)ai for the final state partner produced 
during the backward evolution. The (3 for the on-shell final state partons is 

A'.i = (4.23) 
2aip ■ n 

The p and n vectors are the Sudakov basis vectors for the shower and 



P'±^+l = ZP±^ - PM z'q^ - (1 - z)Ql- (4.24) 



Here is a 2 vector given by (cos0, sin0). By default (p is distributed uniformly in the 
region [0, 27r] but improvements, such as the spin correlations described in chapter Q 

■^If a beam particle was a meson this would just be a (anti-)quark 
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P±i+i = P±i-P±i+i- 



(4.25) 
(4.26) 



All of the ct's and (3's are set until the parton involved in the hard subprocess is 
reached. This fully defines the momentum of the partons. Unfortunately, this does not 
guarantee momentum conservation at the hard subprocess. Instead the momentum of 
the partons must be "shuffled" in order to impose momentum conservation. Rescaling 
the momentum of the partons in the hard subprocess, and correspondingly boosting all 
the partons involved in the evolution of each initial-state partons, allows for momen- 
tum conservation. This rescaling can affect other properties, however, and we want to 
constrain the value of the rescaling so that certain properties are retained. 

In proton-proton collisions we want to conserve the rapidity, Y, and the cm. energy 
squared, M^, of the hard process, while rescaling each incoming parton independently. 
Each parton has its momentum shifted by ki and k2, given in the Sudakov base by 

qi = haiPi + -;-n2 + Pi±, 

012 

q2 = -rP2 + k2p2ni + P2±- (4.27) 
k2 

These shifts, ki and k2, are applied so that the virtuality of the partons is conserved. 
Conserving rapidity and the cm. energy squared requires 

(^ + k2P2 = k,a, + ^, (4.28) 
ki k2 

, , {pi±+P2±y f, , A ^ .PA . . 

Both of these give a quadratic solution for each k value. This leads to four different 
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combinations of values which must be considered. Two combinations of these four have 
negative fc's. This corresponds to flipping the event, e.g. the parton which started 
from the —z direction ends up on the +z direction and vice versa. In some cases the 
reconstruction of the shower fails. In these cases the shower is vetoed and started again 
from the on-shell incoming partons. 



4.3.4 Final-State Shower 

The final-state shower is similar to the initial state shower. The partons start as on- 
shell partons but are taken off shell by the evolution. The initial scale of a parton is 
set according to its shower partner and the partons evolve down to the scale qc in the 
same manner as the initial-state shower, except the probability of emission is modified. 



Instead of Rl^ we have 



RL = \ - (4-30) 

This means that the extra veto, from ()4.1(jj) . is not applied to final-state branchings. 
Though the generation of the branchings is similar, the kinematics of the final-state 
shower is different from that of the initial-state shower. From the results in chapter |21 
we find the relative transverse momentum as 



P±^\ = J{l-Ziyzf{qf-fi')-z.Ql, (4.31 



= ^Jzf{l-z,)^qf-f,^ (4.32) 

for quark branching and gluon branching respectively and /i = max(ma, Qg). 

All of the partons created during the evolution also shower starting at the q that 
they were produced at. Every parton showers until the condition of no more branching 
is met. Once this is reached the parton is put on its mass shell so that momentum 
reconstruction can be done. 
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The partons produced in the hard subprocess have their a set to unity and their (3 
set to 0. Each of children then have their Sudakov variable a set to zai or (1 — z)ai, set 
by the convention of the splitting function. The (3 for the parton corresponding to the 
z momentum fraction is 

Pi+l — 7. 1 {^.66} 

2ai+ip ■ n 

where p and n are the Sudakov basis vectors of the shower. The (3 corresponding to the 
(1 — z) momentum fraction is just 

/3:+i=A-A+i. (4.34) 



Since the final partons are on their mass shells, the initial partons aren't. This means 
instead of having = rn?- they have acquired some virtuality, = g|. The original 
momenta, Pj = {^\Jp] + ^pPj^ the centre-of-mass frame, define some properties of 
the hard subprocess. In the case of e'^e' annihilation we want to preserve the centre-of- 
mass energy 

n 

= E + P] (4-35) 

while keeping the sum of momenta equal to zero. This requires momentum reshuffling. 
To conserve ^/s we rescale the momentum of each jet by a common factor, k, that is 
determined from the equation 

n 

= E \f^+^l (4-36) 

This effectively creates a Lorentz transformation which is applied to all partons in the 
final state. For n = 2 outgoing particles from the hard process ()4.36|) can be solved for 
k explicitly. For n > 2 this is done numerically. 
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We saw earlier that a hard correction is added in to populate the 'dead region' of the 
{n + l)-body matrix element. The parton shower can also be improved in the shower 
regions of phase space by restricting relatively hard gluons that are produced in the 
shower. These gluons are no longer in the domain of of validity of the quasi-coUinear 
approximation. 

Each of the p± values that are generated during the evolution of one jet is tracked. 
If a new p± is larger than any previously generated during the evolution a soft matrix 
element correction is applied to it [73]. This correction assumes all other emissions are 
infinitely soft and we can treat the emission as part of a three-body phase space (e.g. 
qqg). This allows us to compute the three-body variables, (x, x), from the parton shower 
variables (g, z) and the respective Jacobian. The ratio of the hard matrix element and 
the shower is then compared to a random number and the emission is vetoed if the 
ratio is smaller. This requires that the shower approximation is larger than the matrix 
element everywhere in phase space. Otherwise the ratio must have some factor applied 
to ensure that the ratio is always less than 1. This has been studied for several cases in 
chapter |2l and the relevant ratios have been derived. 

4.3.6 Parameterization of Qg 

The cutoff Qg is introduced to regularize the soft gluon singularities in the splitting 
functions. The relative transverse momentum, p± is related to the Sudakov variables of 
the parton branching by 

P± = q±i+i - zq±i. (4.37) 
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Figure 4.2: Available phase space of light (left) and 6-quarks (right) for q ^ qg splitting 
for various values of Qg and depending on the parameterization in terms of S, eq. ()4.4()|1 . 
The dashed lines on the right correspond to a Qg which is not mass dependent, obtained 
by setting m = in ()4.40|) . 

z is required to correspond to a real value of p±. For a gluon splitting this is explicitly 



z-<z<z+, z± = ^ M ± ^1- ^ , , 



(4.38) 



with q > Afi. For quark splittings z is the solution of a cubic but is always in the range 



A* 1 Qg 



(4.39) 



This allows z to be generated within these regions and simply rejected if it lies outside 
phase space. 

Qo is parameterized according to 



6 — 0.3m 
Z3 ' 



(4.40) 



where 6 is the parameter cutoff KinScale in the class ShowerVariables and m is the 
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mass of the parton splitting. This cutoff is used to give the gluons a minimum virtuahty 
which ensures that they are able be put on a mass shell with non-zero mass. As we can 
see from figure 14. 2| the form of this function is also to ensure that b quarks don't have 
an extra artificial cutoff that over-restricts the phase space of splitting b quarks. 



4.4 Hadronization 

Chapter El described in detail the hadronization scheme used in Herwig++. This section 
will discuss the features of the code and the implementation of the specific features. 
Table 14.11 is a table with all the relevant parameters of the hadronization. 

The main driver of the code is the ClusterHadronization class. This class or- 
ganizes and directs which classes are used next. This class also handles the unusual 
case of a hadronization veto. This can occur during the light cluster reshuffling, which 
will be described below. In order to direct the algorithm this class has references to 
PartonSplitter, ClusterFinder, ColourReconnector, ClusterFissioner, and the 
two decayer classes LightClusterDecayer and ClusterDecayer classes. The class 
ColourReconnector by default doesn't do anything. This class is designed to allow 
for a different colour configuration in the forming of the clusters, for example this is 
done in SHERPA [32]. 

The algorithm begins by taking the gluons, which are on a mass-shell, and non- 
perturbatively splitting these into qq pairs. The possible flavours depend on what the 
mass of the gluon is. This is done for all final-state gluons produced in the final-state 
shower (including the final-state evolution of the gluons produced during the initial-state 
shower). The splitting is done by the class PartonSplitter. This class has a reference 
to the GlobalParameters object that defines all the global parameters used throughout 
Herwig++. 
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Parameter 


Class 


Description 


Default Value 


dmax 


ClusterFissioner 


The maximum cluster mass 


3.2 






before a cluster fissions. 






ClusterFissioner 


The cluster mass exponent that 


2.0 




controls cluster fissions. 




PSpltl 


ClusterFissioner 


This is the mass splitting 


1.0 






parameter for udsc flavours. 




PSplt2 


ClusterFissioner 


This is the mass splitting 


1.0 






parameter for b quarks. 




BlLim 


LightClusterDecayer 


Parameter to set the limit of 


0.0 






light b clusters before 








forcing decay to one hadron. 




ClDirl 


ClusterDecayer 


Flag to turn on or off the 


1 






the smearing for non-6 quarks. 




ClDir2 


ClusterDecayer 


Flag to turn on or off the 


1 






the smearing for b quarks. 




ClSmrl 


ClusterDecayer 


Gaussian smearing of non-6 quarks 


0.0 


ClSnir2 


ClusterDecayer 


Gaussian smearing of b quarks 


0.0 


Pwtd 


HadronSelector 


Weight of d quarks 


1.0 


Pwt„ 


HadronSelector 


Weight of u quarks 


1.0 


Pwtg 


HadronSelector 


Weight of s quarks 


1.0 


Pwtc 


HadronSelector 


Weight of c quarks 


1.0 


Pwtfo 


HadronSelector 


Weight of b quarks 


1.0 


Pwt^,^ 


HadronSelector 


Weight of diquarks 


1.0 


SngWt 


HadronSelector 


Weight of baryon singlets 


1.0 


DecWt 


HadronSelector 


Weight of baron decuplets 


1.0 


DKMode 


HadronSelector 


Which hadron decay method to use 


2 



Table 4.1: This is a table of all of the relevant hadronization parameters. Most of the 
parameters are discussed in chapter El and the rest are discussed in this section. 
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Once the gluons have all been split into qq pairs clusters are formed out of the colour 
connected partons. This is done by the ClusterFinder class. This class simply searches 
all coloured partons and creates a colour singlet state out of the colour connected partons. 
There is also a feature to work with partons that are created by colour sources or sinks 
which stem from baryon violating processes. This feature has not been tested though and 
this needs to be done once baryon violating processes are included into Herwig++. After 
the clusters are formed, if a ColourReconnector class was defined to change things, it 
would be called to do so at this stage. 

After the ClusterFinder class has created all of the colour singlet clusters, they 
are passed to the ClusterFissioner class. As was described in chapter |Hl the heavy 
clusters are decayed into lighter clusters. This occurs in three ways. If the mass drawn 
for one of the new clusters in not heavy enough to form the lightest possible pair of 
hadrons, given the set of flavours, the decay is instead forced to C H + C where H 
is a hadron and C is a new cluster. H is the lightest hadron of the flavour of the cluster 
whose mass is too light. The mass generated for the light cluster is changed to equal 
the mass of the lightest hadron. If the mass of the cluster C produced in the decay 
would then violate the phase space bounds it is set to have the largest mass available. 
The next special case is when both clusters are too light to form a hadron pair. This 
decay, C Hi + H2, follows the same procedure as the previous one. Each hadron is 
the lightest hadron of the appropriate flavours and the masses are set to the mass of the 
respective hadron. If this violates the phase space bounds, then the cluster isn't decayed 
via this mechanism. In fact, if this still violates phase space it will be handled by the 
LightClusterDecayer later. The last possible way for a heavy cluster to be decayed is 
directly into two new clusters, C ^ C[ + C'2. In any case, if any of the new clusters is 
still too heavy, it also decays via this algorithm. 

After the heavy clusters have been handled, the clusters which are too light are 
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treated. This is done by the LightClusterDecayer. This takes clusters which are too 
hght to decay into a pair of hadrons and turns them directly into a single hadron. It 
may be possible that a cluster has formed which is even too light to decay directly into 
the lightest hadron. In this case the momentum of this cluster is 'reshuffled' with a 
cluster which is nearby in x-space. The algorithm only 'borrows' momentum from a 
given cluster once during the process. If it is impossible to find enough momentum to 
borrow then the event is vetoed and begun again from the ClusterFissioner stage. In 
the rare case that after vetoing several times there is still not a possible configuration 
to decay all the clusters then the whole event is vetoed. 

The last step of the hadronization is turning the clusters into hadrons. This is done 
by the ClusterDecayer class. This has a reference to a HadronSelector class. The 
HadronSelector class has been implemented with a flag, DKMode, which is used to select 
which of the methods discussed in chapter Elto use. By default this flag is set to 2, which 
is the new method. If it was changed to the implementation of the HERWIG method 
would be used. A value of 1 is for the Kupco method. Once a pair of hadrons has been 
selected by one of the methods then the cluster is decayed. The decay products are 
generated in the direction of the constituent quark. This direction can have a Gaussian 
smearing applied to it to. If the flag ClDir is set and ClSmr is larger then 0.0001, then 
cluster decay is in the direction 



cos 6 = cos 6g cos 6. 



smr 



— sin 6q sin 6. 



smr ) 



(4.41) 



where 9q is the 9 of the original constituent quark and 



cos 6. 



smr 



1 + ClSmrlog7^, 



(4.42) 



with TZ a random number. The azimuthal angle is distributed uniformly in [— tt, tt]. 
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In chapter Q the general scheme of the decays was presented. Since the algorithms 
implemented in Herwig++ are very basic, in this section I will explain the way to setup 
the decay modes in Herwig++ and how one would implement a new set of decay modes. 

In Herwig++ the default decays are defined in the file Hw64Decays . in. In this file 
all of the decay modes are included. An entry of a decay mode takes the form 

decaymode a->b, c,d,...; BR 0/1 Class 

for the decay a ^ b + c + d+ . . . , BR is the branching ratio and the 0/1 option specifies 
whether to turn this mode off or on, respectively. The last argument specifies which 
class to use to perform the decay. By default Herwig++ has only the HwDecayer class, 
the HeavyDecayer class and the QuarkoniumDecayer class. All of these classes use a 
parameter MECode to specify the matrix element code to use. 

The class HwDecayer is the most common class used in the default Herwig++ decays. 
This just implements the decays as described in chapter ^ based on the MECode param- 
eter. A value of 100 uses the free massive {V — A){V — A) matrix element, for example 

—>■ UeCUr- 101 is the code to use the bound massive matrix element. This would 
be used for a decay such as K~ — > Uge^Tr'^. Any other value of the parameter uses an 
isotropic n-body decay. 

The class HeavyDecayer takes heavy mesons and decays the heavy parton weakly. 
The W produced in the decay is also decayed all in one step. This produces 4 partons, 
two from the decay of the W, one from the production of the W and the remaining 
spectator parton. For example a chain could look like, B^ dcW^ dcdu. The 
decay products have the colour connections set properly and the administrative class 
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HwDecayHandler directs the program to reshower and rehadronize these partons. 

The last decayer, the QuarkoniumDecayer, is designed to work with heavy mesons 
and baryons that decay into gluons and quarks. Possible modes are to 3 gluon states, 
2 gluons and a photon, just 2 gluons or a qq state. This class handles a MECode of 130 
by using a positronium matrix element. Decays from this class also lead to showering 
and hadronization of the partonic decay products. An example of this type of decay is 
Vc 99- 

Implementation of a new decay mode requires the definition of only two functions. 
A new decay matrix element must inherit the Decayer class from ThePEG. The two 
functions that must be implemented are 

• bool accept (const DecayMode &) const; this function is used to indicate if a 
particular decay mode can be handled. For example if you implemented a decay 
matrix element for decays, if the incoming particle was not a r^, this function 
would simply return false. 

• ParticleVector decay (const DecayMode &, const Particle &) const; this 
function returns the decay products with their momentum set. The input is the 
decay mode selected based on the branching ratios, and the instance of the particle 
which is decaying. 

Improving the decay modes of Herwig++ is currently underway. This is a lengthy 
process and as long as experimentalists study the decays of particles, improving the 
decay mode matrix elements, as well as the branching ratios, will always be possible. 
But as will be shown in the next chapter, even with the simple decay matrix elements, 
fits to particle spectra are good. 



Chapter 5 



Results 



5.1 Introduction 



This chapter presents results from the Herwig++ . The results given here are for e^e~ 
annihilation events, as this is the first step in the redevelopment of HERWIG. In order 
to have full control of the basic physics steps that are simulated, it was thought to be 
very important to put the new generator on a firm basis with respect to LEP and SLC 
results before upgrading it to be able to deal with initial-state showers and the other 
requirements for the simulation of lepton-hadron and hadron-hadron collisions. There- 
fore, thorough tests have been performed on the predictions of the generator against 
a wide range of observables that have been measured at LEP and SLC. We have also 
explored the sensitivity to the most important parameters and cutoffs. These results are 
not the result of a high-precision tuning: the main aim here is rather to show the results 
of the program and that it is able to give results as acceptable as those generated by its 
predecessor HERWIG for a reasonable choice of parameters. 
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As was discussed in detail in the last three chapters, the main stages of the simulation 
of events is the same as in HERWIG [62]. However, in comparison to its predecessor, 
Herwig++ features a new parton shower and an improved cluster hadronization model. 
At present, hadronic decays are implemented in the same fashion as they were in HER- 
WIG. 

As discussed in chapter |31 the program is based on the Toolkit for High Energy 
Physics Event Generation (ThePEG) [64] and the Class Library for High Energy Physics 
(CLHEP) [65]. They are utilized in order to take advantage of the extended general func- 
tionality they can provide. The usage of ThePEG unifies the event generation framework 
with that of Pythia7. This will provide benefits for the user, as the user interface, event 
storage etc. will appear to be the same. The implementations of the physics models, 
however, are completely different and independent from each other. 

The simulation of e~^e~ annihilation events starts with an initial hard process e'^e" 
7*/Z° ^ qq + 77. The final state photons simulate QED radiation from the initial 
state, so that a radiative return can be properly simulated. For these results we are 
only interested in the details of the QCD parton shower in the final state. The final- 
state parton shower starts with a quark and antiquark that carry momenta Pq and Pg, 
respectively, and have an invariant mass squared of = {pq + PqY- For the e~^e~ 
results, the only detail we are concerned with in relation to initial-state radiation is that 
the centre-of-mass frame of the gg-pair is slightly boosted with respect to the collider 
laboratory frame and that Q may be different from the e^e~ centre-of-mass energy. We 
have made sure that the applied cuts on the energy of the annihilating e~^e~ subsystem 
are the same as those used in the experimental analyses. 

Currently, proton-proton collisions are being studied. These results are simulated by 
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starting with the Drell-Yan hard process qq — > 7*/^" £~^i~ , with £^ a charged lepton. 
These results concentrate on the effects of the initial-state radiation from the incoming 
quarks and gluons. 

5.1.1.1 Parton shower 

As has been shown in the preceding chapters, the partonic evolution from the large scale 
of the hard collision process down to hadronic scales via the coherent emission of partons, 
mainly gluons, is simulated on the basis of the Sudakov form factor. Starting from 
the hard process scale Qq, subsequent emissions at scales Qi and momentum fractions 
Zi are randomly generated as a Markov chain on the basis of the soft and coUinear 
approximation to partonic matrix elements. Chapter El has shown that for Herwig++ 
we have chosen a new framework of variables, generically called (g, z). Here, g is a scale 
that appears naturally in the quasi-coUinear approximation of massive partonic matrix 
elements and generalizes the evolution variable of HERWIG to the evolution of massive 
quarks, z is a relative momentum fraction; the evolution is carried out in terms of the 
Sudakov decomposition of momenta in the frame where the respective colour partners 
are back-to-back. As in HERWIG, the use of the new variables allows for an inherent 
angular ordering of the parton cascade, which simulates coherence effects in soft gluon 
emission. The details of this underlying formalism have been described in chapter |21 

The most important parameter of the parton shower that we will be concerned with 
in this chapter is the cutoff parameter Qg, which regularizes the soft gluon singularity 
in the splitting functions and determines the termination of the parton shower. This is 
set by 6 in ()4.40|) . Less important but relevant in extreme cases is the treatment of the 
strong coupling constant at low scales. We have parametrized as{Q) below a small scale 
Qmin > ^QCD in different ways. We keep Qmin generally to be of the order of 1 GeV, 
where we expect non-perturbative effects to become relevant. Below that scale as{Q) 
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• set to zero, as{Q < Qmin) = 0, 

• frozen, as{Q < Qmin) 

• linearly interpolated in Q, between and as(Qmin), 

• quadratically interpolated in Q, between and asiQmin) ■ 

We put the final partons of the shower evolution on their constituent mass shells, since 
the non-perturbative cluster hadronization will take over at this scale, so we usually 
have kinematical constraints that keep Q above Qmim in which case the treatment below 
Qmin is irrelevant. Typically, as{Qmm) ~ 1 here. 

5.1.1.2 Hadronization and decay 

As discussed in chapter |3J the partonic final state is turned into a hadronic final state 
within the framework of the cluster hadronization model of HERWIG [74]. All three 
methods of cluster decays have been implemented in Herwig++, but the new cluster 
hadronization model is used for the results given in this chapter. The emerging hadrons 
are possibly unstable and eventually decay. The decay matrix elements and modes 
correspond to those in HERWIG. 

5.2 e+e~ Annihilation 

This section presents the results for e+e~ annihilation events. The properties of differ- 
ent measurements are discussed and the comparison of Herwig++ to data is presented. 
Histograms for all the distributions have been booked in the same bins as the experi- 
mental data. For a given bin i we then compare the data Di value with the Herwig++ 
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Parameter 


Herwig++ 


Initial 




0.118 


0.114 


S/GeV 


2.3 





mJGeV 


0.750 





<5min/GeV in ^^(Qmin) 


0.631 





Cl„ax/GeV 


3.2 


3.35 


CXpQ^ 


2.0 





PSpltl 


1 





PSplt2 


0.33 


— 


BlLim 


0.0 




ClDirl 


1 




ClDir2 


1 




ClSmrl 


0.40 




ClSmr2 


0.0 




Pwtrf 


1.0 




Pwt„ 


1.0 




Pwts 


0.85 


1.0 


Pwtc 


1.0 




Pwtft 


1.0 




Pwtgg 


0.55 


1.0 


Singlet Weight 


1.0 




Decuplet Weight 


0.7 


1.0 



Table 5.1: The parameters for Herwig++ used in this study. The first group are shower 
parameters, the second are all of the hadronization parameters. In the third column we 
show initial values of our study, taken from HERWIG. 
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Monte Carlo result Mj. Given the data errors SDi (statistical plus systematic, added 
in quadrature) and Monte Carlo errors SMi (statistical only), we can calculate a ^or 
each observable. We keep the statistical error of the Monte Carlo generally smaller than 
the experimental error. In distributions where the normalization is not fixed, such as 
momentum spectra, we allow the normalization of the Monte Carlo to be free to min- 
imize x^. The normalization is then tested separately against the average multiplicity. 
In all other cases we normalize histograms to unity. 

As we do not want to put too much emphasis on a single observable or a particular 
region in phase space where the data are very precise, in computing we set the relative 
experimental error in each bin to max{6Di/Di,5%). This takes into account the fact 
that the Monte Carlo is only an approximation to QCD and agreement with the data 
within 5% would be entirely satisfactory. The general trend for the preferred range of a 
single parameter was however never altered by this procedure. 

After normalization the ratio 



is computed for each bin in order to see precisely where the model fails. This ratio as 
well as the relative experimental error and the relative contribution of each bin to the 
of an observable is plotted below each histogram. 



We have taken x^ values for hadron multiplicities into account in the same way as we 
weighted the event shapes. In general the multiplicities of individual particle species are 
sensitive to a completely different set of parameters. The general strategy was to get a 
good value for the total number of charged particles with a reasonable set of values for 




(5.1) 



5.2.1 Strategy 
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the parton shower cutoff parameter S and the maximum cluster mass parameter Clmax- 
Once this was fixed, the hadronization parameters that determine the multiphcities of 
individual particle species were determined. Following this we compared this 'preferred' 
set of parameters with the 'default' set from HERWIG. The resulting parameter set is 
shown in Table ISTTl 

A wide range of observables have been tested in order to study the aspects of the 
model. Event shape variables and multiplicities are considered in order to test the 
dynamical aspects of the parton shower and hadronization models, which are closely 
linked at their interface by the parton shower cutoff parameter, 6. Ideally, the models 
should combine smoothly at scales where ~ 1 GeV. Many of the figures shown in 
this chapter contain three sets of plots per figure. In order from top to bottom these are 

• the actual distribution. The Herwig++ result is plotted as a histogram together 
with the experimental data points; 

• the ratio Ri together with an error band showing the relative statistical and 
systematic errors; 

• the relative contribution of each data point to the total of each plot. 



5.2.2 Hadron multiplicities 

The charged particle multiplicity distribution and the overall multiplicities of a wide 
range of hadron species have been taken to test the overall flow of quantum numbers 
through the different stages of the simulation. This also allows a thorough test of the 
new hadronization model, developed in chapter El against the measured observables. 

Table ini21 showed the results of the new cluster hadronization algorithm in comparison 
to the old algorithm. Even before systematic tuning, we can see that the overall results 
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1 I 
Herwig++ 1.0 

6 = 1.7 GcV 

S = 2.3 GeV 

S = 3.2 GeV 

OPAL 99 




Figure 5.1: The distribution of the charged particle multiphcity. 

are in better agreement with the data than those of HERWIG, with fewer results that 
differ from the data by more than three standard deviations (indicated by a star in 
the table). We can also see from table ESI that the difference between the models is 
quantified by their x^- 

A very well measured property, and therefore important to get accurate, is the dis- 
tribution of charged particle multiplicity. Figure 15.11 shows the results of Herwig++ 
compared to OPAL data [75] and is found to be in fairly good agreement. There is 
an excess of lower multiplicity events, however. It is also shown that varying 6 doesn't 
greatly alter this distribution, another confirmation that the interface between the par- 
ton shower and the hadronization is consistent. 
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This measurement is the multiphcity of (mini-)jets in e"'"e~-colhsions for different values 
of the jet resolution ?/cut- We use the Durham- or /c_l -clustering scheme [76] throughout 
this chapter for jet observables. To be specific, for a given final state the jet measure 

2mm(Ef,E^) 

is calculated for every particle pair (^, j). The particles with minimal 'distance' are 
clustered such that the momentum of the clustered pseudo-particle is the sum of the 
four-momenta of the constituents. The jet multiplicity is then the number of pseudo- 
particles remaining when all yij > ?/cut- This inclusive observable has been predicted and 
measured at LEP energies and will test the dynamics of the parton shower as well as the 
interface between parton shower and hadronization. We use the Kt Jet-package [77] that 
implements the above jet-finding algorithm in C++ and have written a simple wrapper 
around it in order to use it with the particle record of Herwig++. 

Figure shows the average number of jets (njets) at the Z^-pole, as a function of 
the Durham jet resolution, i/cut, for various values of the cutoff parameter 6. At the 
parton level (top left) the jet multiplicity varies substantially toward smaller values of 
i/cnt, saturating at the number of partons that are present in a single event. The order of 
magnitude of the visible saturation scales is characterized for each flavour by different 
cutoff values Qg as ysat = Q^/Q^- For example, at Q = 91.2 GeV and 6 = 2.3 GeV, the 
saturation scale for light quarks is of the order 10~^ while for 6-quarks it is of the order 
10"^. 

During hadronization, low parton multiplicities lead to large mass clusters which, 
as described before, tend to decay into low mass clusters below the cutoff mass, Clmax- 
This has been fixed to its default value for the results given in this chapter. Figure 15.21 
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(top right) shows that the hadronization compensates for lower partonic multiphcities, 
giving a result which is insensitive to 6 at the hadron level. This means that we have a 
smooth interface between the perturbative and non-perturbative dynamics of the lower 
end of the parton shower and the cluster hadronization model. On the hadron level we 
describe LEP data from OPAL [78] well. 

In order to test the sensitivity of Herwig++ against the variation of the centre-of- 
mass energy we calculate the jet multiplicities at PETRA and LEP II energies as well 
as LEP I energies ()5.2[ bottom). The comparison to JADE [79] and OPAL [78] data 
shows a good agreement. For the generation of all the Monte Carlo data we applied the 
same cutoffs on the energy of the partonic subsystem as was done in the experiments. 

The other curves in figure 15.21 show the prediction for the jet multiplicity [80] from 
the resummation of leading logarithms. Note that the parameter Aqcd in the resummed 
calculation is not Ay^. For a value of Aqcd = 500 MeV we can see that there is good 
agreement with the data and the Herwig++ results throughout the perturbative region, 

Vent > lO-\ 



5.2.4 Jet fractions and 

These measures give a closer look 'into' the jets. This is done by considering the rates 
of jets at a given value of jjcut in the Durham scheme. The jet fraction is given by 

Rn = (5.3) 

forn = 2 up to n = 6 jets. Presented here are also the distributions of 1^, the ?/cut-values 
at which an n+l-jet event is merged into an n-jet event in the Durham clustering scheme. 
The results are presented here for n = 2 up to n = 6, without n = 5. These distributions 
will not only probe the dynamics of the parton shower but also the hadronization model; 
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Figure 5.2: Jet multiplicities for different values of the cutoff parameter 6 and different 
centre-of-mass energies. 



CHAPTER 5. RESULTS 



149 



at the lowest values of ?/cut ~ {qc/QY the dynamics is dominated by the hadronization. 

Figure ESI compares the results from Herwig++ with LEP data from [78] and shows 
good agreement. On the hadron level these predictions are not very sensitive to the 
cutoff parameter S. The results of the are not shown here but show similarly good 
agreement. The -Re plot contains all jets n = 6 and greater. 

The Durham y„ distributions are given in figure These are histograms of the ?/cut 
values at which the n + 1-jet event in the Durham jet clustering scheme is merged into 
an ra-jet event. This resolved more of the internal structure of the jets than the n-jet 
rates alone. Overall, the agreement between the model and the data is good. There is 
a tendency to exceed the data at low y„. This is a problem that was also present in 
HERWIG. 

5.2.5 Event shapes 

Event shape distributions have been measured to very high accuracy at LEP and aim at 
resolving the properties of the parton shower quite thoroughly. All of the event shapes 
given here, except 5*, P and A, defined below, are 'infrared safe'. This means that they 
can be computed in perturbation theory. In order to test the dynamics of the parton 
shower in Herwig++ in more detail we consider a set of commonly used event shape 
variables. Not only the collinear region of the parton shower is probed in greater detail 
but also the regions of phase space which are vetoed as matrix element corrections. We 
compare all results to DELPHI data [81]. 

The thrust is a well studied property. The definition of the thrust is given by 



T = max 



n 



(5.4) 
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Figure 5.3: Jet rates in the Durham algorithm for different values of the cutoff 6 
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Figure 5.4: Durham Yn distributions for different values of the cutoff S. 
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Finding the vector n is a computational intensive task. This can instead be reduced 
using physical arguments to a simpler procedure. We start by separating the sum into 
two different parts: those where Pa ■ n > and those where ■ n < 0. If we first take 
out the normalization factor 

N = J2\Pal (5.5) 

a 

we get 

T = ^ max I ^ \pa-n\+ ^ |p„ • n| j . (5.6) 

\PaTl>0 Pa-n<0 / 

The magnitudes of the dot products can be removed and this gives 



T = — max ( 



Pa-n- ^ Pa-n] . (5.7) 

Tl>0 Pc,Tl<0 / 

Since the n is independent of the sum it can be taken outside the summation yielding 



T = ^max(P+(n) - P_(n)) ■ n, (5. 

N n 



where P^{n) is the sum of all momenta in the same hemisphere as n and P^{n) is the 
sum of all momenta in the other hemisphere. Momentum conservation says (in the cm. 
frame) that P+ = — P_ so the thrust is given by 

2 

T = — max P+(n) ■ n. (5.9) 

From this it is obvious that the maximum value for a given n is when the vector 
lies parallel to P+(n). Previously, to find the vector n that maximizes ()5.4j) we would 
have to consider all possible combinations of momenta, which has complexity of order 
A^!. After deriving ()5.9|1 this has been reduced to considering all sets of two vectors 
which define a plane. This plane divides the space into two hemispheres from which the 
thrust can be computed by simply summing all the momenta in one of the hemispheres 
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Figure 5.5: Thrust without (top left) and with (top right) matrix element corrections 
switched on, thrust major and thrust minor (bottom). 
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and taking the magnitude. Since the two vectors chosen to define the plane do not 
unambiguously lie in either hemisphere, each possibility must be considered. This has 
reduced the problem to simply iterating over all sets of two momenta and recording 
which set produces the maximum. This is now only order A^^. 

There are two other measures of the data related to the thrust. One is called the 
thrust major and the other the thrust minor. These are defined as 



M = max , , ' , (5.10) 

m = — ^^—j — ■ — , (5.11) 



where M is the thrust major, m is the thrust minor and = rix x um- 

The thrust is a measure used to describe how 'pencil-like' the event is. High thrust 
means that the event is more 2-jet like, where as lower thrust means that the event is 
much more planar or spherical, thus it has more than 2 jets. The thrust major is used to 
describe the major component of the momentum in a plane perpendicular to the thrust 
axis. High values of the thrust major are usually indicative of planar events. The thrust 
minor, therefore, is used to describe the remaining degree of freedom of the momenta. 
High values of the thrust minor correspond to spherical events. These are events that 
have at least 4 jets. 

In fig. 15.51 we show the distribution of thrust and thrust-major and thrust-minor. 
These variables are all obtained from the equations given above. The thrust distribution 
is shown with and without matrix element corrections switched on. The prediction 
without matrix element corrections is very much better than that of HERWIG, owing 
to the improved shower algorithm. It is interesting that the matrix element corrections 
seem to generate almost too much transverse structure, leading to event shapes that are 
less two-jet-like. On the other hand, there is also a slight excess of events close to the 
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Figure 5.6: C parameter and D parameter distribution. 



two-jet limit. 



Though the thrust describe 2-jet hke events quite well, multi-jet events are less un- 
derstood in terms of these measurements. The C and D parameters are used to describe 
three- and four-jet-like events. These are given by combinations of the eigenvalues of 
the linear momentum tensor 



Ea \Po\ 



(5.12) 



The definition of C and D is then 



C — 3 (A1A2 + A2A3 + A3A1) , 
B = 27A1A2A3, 



(5.13) 
(5.14) 
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Ai + A2 + A3 = 1. 



(5.15) 



Both of these parameters have a coefficient defined so that the range of the parameter 



It is remarkable how well distributions like C and D parameter (fig. 15. 6|) which are 
sensitive to three- and four-jet-like events are described by our model even though we 
are limited to three jet matrix elements plus showers. Here again we have in fact a small 
excess at high values. 

We show also in fig. 15.71 the distributions which are obtained from a quadratic mo- 
mentum tensor 



The three measures that arise from this are sphericity (S) , planarity (P) and aplanarity 
{A). These are given by 



is [0,1]. 





(5.16) 



y 



s 



I (A2 + A3) , 



(5.17) 



P 




A 




(5.19) 



where Aj is the ith eigenvalue and they obey the relations 



Ai > A2 > A; 



Ai + A2 + A3 = 1. 



(5.20) 



These distributions put more emphasis on high momenta. As the names imply these 
distributions indicate the events that are spherical, planar or aplanar. The sphericity 
axis is just = Vi where Vi is the eigenvector of the ith eigenvalue. This is made 
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Figure 5.7: Sphericity, planarity, and aplanarity parameter distribution. 
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Figure 5.8: The wide and narrow jet broadening measures -Bmax and i?n 



obvious by seeing that sphericity is really defined as 



Ai). Therefore, events which 



have high sphericity have momentum which tend to be along the sphericity axis, n^. As 
was the case for the thrust-related distributions, we tend to have slightly wide events. 

In addition we consider the jet broadening measures -Bmax and -Bmin and the hemi- 
sphere jet masses (fig. 15.81 and fig. 15. 9^ . The jet broadening measures are defined 
as [82,83] 

(5.21) 



5rnax = max ' I — . 



where ut is the thrust axis and Hi indicates one of the two hemispheres defined by the 
plane normal to ut- If Pk ■ f^T > then pk is in hemisphere. Hi, otherwise it is in 
hemisphere H2. -Bmin is then 
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Figure 5.9: The high and low hemisphere masses. 



-Bmax, -Bmin measure the scalar momentum transverse to the thrust axis for the wider 
and narrower jet hemispheres respectively. We can see from figure 15.81 that there is 
good agreement between the model and the data. We have also looked at two other 
jet broadening measures, Sdis = -Bmax — -Bmin and Sgum = -B^ax + B^^^^. These are not 
shown here as they contain the same information as i?max and i^min- 

The last jet measure we show here is that of the hemisphere masses. These also use 
the same definitions of the two hemispheres with respect to the thrust the jet 

broadening measures, but these measure the total momentum squared within a hemi- 
sphere. The high hemisphere mass is 



high 



1 max ( Y: Pk] 



(5.23) 
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and the lower hemisphere mass is 



- mm i y^Pk] 

S i=l,2 \ ^ ] 



Miow = T mill I 2^ Pk] , (5.24) 

where s is the cm. energy squared. In both cases we can see from figure EUl that the 
agreement between model and data is good. 

5.2.6 Four jet angles 

We show the four-jet angles in fig. 15.101 They are considered only for events where we 
have a four-jet event at ?/cut = 0.008. Each of the different angles measures a property 
of the four-jet configuration. 

If we denote the cosine of the angle between two vectors, a and b as C{a, b), we can 
define the four different jet angles as 



cosxbz = C ((pi X Pa), (P3 X p^)) , (5.25) 

«! + a2 
_ 2 _ 

cosGnr = C((pi -p2),(P3-P4)), (5.27) 



cos $KSW = COS 



(5.26) 



cosa34 = C(p3,p4), (5.2^ 



where 



cosai = C ((pi X P4), (p2 X P3)) , (5.29) 
cosa2 = C ((pi X P3), (p2 X P4)) . (5.30) 



Pi is the three-momentum of the hardest jet and p4 is the three-momentum of the softest 
jet. 
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Despite the fact that we do not have any matching to higher order matrix elements, 
as was proposed in [27] and implemented in [84], the agreement between model and 
data [85] is remarkably good. We expected the implementation of hard and soft matrix 
element corrections in Herwig++ to improve the description of these observables but 
we did not find very significant differences with or without the application of matrix 
element corrections. 



5.2.7 Single particle distributions 

Using the thrust axis defined in ()5.4|) . one can define the transverse momentum of a 
particle with respect to this axis. We can then define a plane given by the thrust axis and 
the thrust major axis, um as the event plane. Therefore, the plane given by the thrust 
and the thrust minor, n^, is considered out of the event plane. In fig. l5.11l we have single 
charged particle distributions. The distributions are of the transverse momentum within 
the event plane, j^^, and the transverse momentum out of the event plane, The 
momenta in the event plane are shown with and without matrix element corrections. In 
contrast to the thrust distribution we find that the matrix element corrections actually 
improve the distribution. Furthermore, and the rapidity along the thrust axis, 

i/T, are rather well described. A similar technique can be used to define the transverse 
momentum with respect to the sphericity axis n^. We do not show these but they have 
similar features to the distributions with respect to the thrust axis. 

As we saw in chapter IHl we can consider the distribution of scaled momentum Xp = 
2\p\/Q of charged particles. The results in chapter El were given to show that the 
hadronization technique could describe these distributions quite well. In fig. 15.121 these 
distributions are given to show their dependence on the shower cutoff, 6. In addition to 
the full distribution we also consider the results from light {uds), c and b events. ^ In 

''The flavour of tlie quark-antiquark produced in the initial hard process 
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Figure 5.11: Momentum distributions of charged particles with respect to the thrust 
axis, (with and without matrix element corrections), and y'^. 
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all cases we compare with data from SLD [86] . The charged particle distribution is well 
described in all four cases, in fact somewhat better for heavy primary quarks. 



5.2.8 Identified hadron spectra 

As in the case of all charged particles we can compare identified particle spectra from 
events of different flavour to SLD data [86]. Data for tt^ (not shown, being almost 
equivalent to all charged particles), and {p,p) is available. In fig. 15.131 we see the 
data for {p, p) spectra from events of different flavour. For large values of Xp we clearly 
overshoot the data in light flavoured events, producing the 'bump'. This is somewhat 
compensated by the heavy quark events which in turn seem to prefer lower values of Xp. 
The origin of the 'bump' is not well understood. We believe that this feature is related 
to the hadronization, being similar to but smaller than that seen in HERWIG, but this 
has not been shown conclusively. 

Fig. 15. 141 shows distributions for and A, A. Both are rather better described than 
the proton spectra but the distribution of A, A tends to have a similar, though smaller, 
'bump' in comparison to data from ALEPH [87]. 



5.2.9 B fragmentation function 

The hadron fragmentation function is defined as the distribution of 



xe = (5.31) 



-'cm 



where Eh is the energy of the hadron. Using this we can look at the fragmentation 
function of different flavours or species of hadrons. In fig. 15. 151 we consider the B hadron 
fragmentation function in comparison to data from SLD [88]. This is the fragmentation 
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Figure 5.12: The scaled momentum distribution Xp of charged particles for all events as 
well as for uds, c and b events separately. 
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Figure 5.13: The scaled momentum distribution Xp of protons, shown separately for all 
events as well as for uds, c and b events. 
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Figure 5.14: Distribution of scaled Kaon momentum and A, A momentum. 

function for all weakly decaying B hadrons. We have also considered data from ALEPH 
[89] for comparison, though these are not shown here. We can describe the data quite 
well without any additional tuning of the hadronization model to this data. The parton 
shower formulation in terms of the new variables [36] and taking quark masses in the 
splitting functions into account clearly improves the description of heavy quark events. 



From figure IHTTKl we tend to bias the fit towards the b of 2.3 GeV, as the improved 
treatment of the h quarks was the main motivation for deriving the new variables and 
using the massive splitting functions. HERWIG couldn't describe the data as well 
even with the extra flags added to the hadronization model to parameterize B hadrons 
differently. 
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Figure 5.15: The i?-hadron fragmentation function. For different values of the cutoff 5. 
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In Tab. 15.21 we show a list of values for all observables that were studied during our 
analysis, including those not shown in the plots. The most sensitive parameters were 
the cutoff value 6 and the use of (hard plus soft) matrix element corrections. The table 
shows three values of 6: our preferred value of 5 = 2.3 GeV as well as the lowest and 
highest values that we considered. 

The results should be interpreted with care. The overall trend suggests that we should 
prefer a large cutoff scale. However, we have just averaged over all possible observables. 
Taking a closer look, we may want to weight different observables in a different way. 

In more detail, the general trend is the following: event shapes, jet rates and differ- 
ential jet rates prefer a low cutoff. The single particle distributions along the thrust and 
sphericity axes prefer a small cutoff value. The ynm distributions prefer either a high 
or a low cutoff value. The spectra of identified particles tend to prefer the high cutoff 
value with some exceptions for light quark events. The B fragmentation function clearly 
prefers the intermediate value. 

In addition, as indicated in sect. 15.2.1} we found that the measured yields of identified 
particles clearly prefer the value 6 = 2.3 GeV. 

5.3 Conclusions 

We have achieved a complete event generator for e+e^ annihilation into hadrons. The 
main physics features, in comparison to the previous versions of HERWIG, are an im- 
proved parton shower, capable of properly describing the perturbative splitting of heavy 
quarks, and an improved cluster hadronization model. 

We have tested our model against a wide range of data from e+e" colliders and are 
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Table 5.2: values for all observables we studied and a relevant subset of parameters. 
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able to give a good general description of the data. 

For many observables the description of the data has been improved with respect to 
HERWIG. The new parton shower has a number of remarkable features. The need for 
matrix element corrections has decreased. The main reason for this is the use of improved 
splitting functions, which give a far better approximation of the matrix elements in the 
region of collinear gluon emissions. We can describe observables involving light or heavy 
quark splitting with a unique set of parameters. The new hadronization model also 
improves the description of identified particle spectra and multiplicities. 

The detailed analysis of our results leaves us with a recommendation: the set of 
parameters that is shown in table 15. Il This set of parameters is understood as a weighted 
compromise in order give a good overall description of the data we have considered so 
far. We did not aim at a complete tuning of the model, but rather wanted to study its 
ability to describe the broad features of the data, which turned out to be very successful. 

Work is currently under way testing the parton shower on initial state radiation. A 
model for the soft underlying event in hadron-hadron collisions is also under develop- 
ment. The aim is to have the code tested and debugged so a complete event generator 
for the simulation of Tevatron and LHC events is available. 



Chapter 6 

Effective Potential Analysis: 
Effective 

6.1 Super symmetry 

Though the SM presently describes physical phenomena extremely well there is reason 
to believe that it will eventually be insufficient to describe physics at higher scales. The 
scale tested to date is of the order of 100 GeV. The Planck scale, Mp = (SvrG'Ncwton)^^^^ = 
2.4 X 10^*^ GeV, is the scale at which gravitational effects will be of the same order as 
the other interactions. This scale is 16 orders of magnitude higher than the currently 
observed phenomena. On an intuitive basis this is suggestive that there must be new 
phenomena occurring between these scales. This is commonly referred to as the hierarchy 
problem. 

There is also another problem, known as the fine tuning problem. This is due to the 
Higgs boson. The mass of the Higgs boson is restricted to be of the order 100 GeV. 
Figure IFTD shows the diagrams that correct this mass at one loop. From fig. Ifi.lk we get 
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/ 

(a) (b) 

Figure 6.1: These diagrams are the one loop corrections to the Higgs mass for fermions 
and scalars. 

the correction due to Dirac fermions with mass rrif [90] 

^"^H = tA [-2A?/y + SmJ In {Auv/mj) + ...], (6.1) 

where Auv is an ultraviolet momentum cutoff at the scale of new physics and A/ is 
the coupling derived from the term — A//$/. This problem directly affects only the 
correction to the Higgs scalar boson (mass)^ because the fermions and gauge bosons 
do not have the quadratic sensitivity to Auy. If Auv is of the order of Mp then this 
correction is about 30 orders of magnitude larger than the expected mass-squared of 
the Higgs boson. This will indirectly affect the SM particles and their known masses 
as corrections to these masses have some dependence on mn- Therefore, this cannot be 
the case. 

The Higgs mass also receives a contribution from fig. 16.1b from the scalar particles. 
For a real scalar field, this contribution takes the form [90] 

Am^ = [Alv - 2"^| In (Auv/ms) + ...]. (6.2) 

Due to the opposite signs in ()6.2|1 versus 1)6.11) there could be a particular combination 
of particles and masses that exactly cancel the divergences. Though plausible, it seems 



H 



H 
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extremely unlikely that all of the parameters that are currently in the SM plus any 
particles occurring with a higher mass are tuned exactly to cancel this. 

Instead a systematic cancellation of the divergences can occur. We see that if we 
have two complex scalar fields for each fermion field which couple to the Higgs boson 
exactly as Xs — |A/| then they exactly cancel the quadratic terms. This symmetry 
between fermions and bosons is known as supersymmetry [90-92]. 

Of course this symmetry must be broken because we do not observe bosonic states 
of the same mass as a fermion partner. The scale at which it is broken is Msusy and 
the SUSY particles will have masses of the order of Msusy- In SUSY models [90] this 
value is expected to be at most 1 TeV in order to allow for a Higgs VEV which will give 
the correct values of mz and rnw- 

At a scale of 1 TeV, the experiments at LHC should be able to discover SUSY 

particles. Research is still ongoing into new ways to use the general purpose experiments 
ATLAS [93] and CMS [94] to search for SUSY particles (for example [95,96]). 

6.1.1 Superpotential 

The fermion and boson that are supersymmetric partners are combined into a multiplet 
called a supcrmultiplct. There arc two kinds of supermultiplcts, a chiral and a gauge 
supermultiplet. A chiral supermultiplet is a supermultiplet formed by the matter fields. 
The gauge fields and their superpartners form gauge supermultiplcts. These are the 
objects that enter into the supersymmetric Lagrangian. This in turn can be written in 
terms of the superpotential, which will be described in this section. 

It can be shown that the most general set of renormalizable SUSY interactions for a 
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field ipi and its SUSY partner can be written as 



Ant = --W''i^^iJJ + W'F, + C.C., (6.3) 



where the Ws are explained below, the F is an auxiliary term and the c.c. indicates 
complex conjugate. The auxiliary terms are used as a bookkeeping device and are 
eliminated by the equations of motion. Doing so shows that Fi = —W*. 

We can now define the superpotential as 

W = ^M^'M, + ^y'^^Mk, (6.4) 

where M*-^ is a mass matrix for the fermions and y'-''^ is a Yukawa coupling of scalar, 0^, 
and fermions ipiy'^j- This is related to the W's in eq. (j6.3|) by 

6(f)i6(f)j 

W, = (6.6) 



By replacing the F terms and adding the kinetic term, we find the complete La- 
grangian for a chiral supermultiplet is 

-Cchiral = Cun " ^ (w'^^.^j + VT^^I^j) - V{<j), <f)*), (6.7) 

where V{(j), </>*) = W^W* is the scalar potential and £kin is the standard kinetic terms 
for a fermion and a scalar. 



The gauge field interactions also have to be written in terms of their superpartners. 
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For a gauge boson, A" and its superpartner fermion, A°, this is 



where is the same term as in the SM and D"" is another auxihary field, cr^ and cr'^ 
are given by 

a'^ = (l,aO, a^ = (l,-a^), (6.9) 
and are related to the standard 7-matrices by 



7" 



^0 ..^ 



(6.10) 



It can be shown that there are two other terms that can be added and the La- 
grangian will remain invariant under supersymmetric transformations. When added to 
the Lagrangian they give 

C = C,ur.i + -Cgauge " V2g [(0*T» A'^ + A^^^ (^^TV)] + g (0*T» D^ (6.11) 

where T° is the generator of the group A" is a mediator of. It is from this equation 
that we can find the equations of motion to remove D"^. We find that D"^ = —g {(f)*T°'(f)) 
and the ^D'^D"^ from ()6.8p combines with the last term in ()6.11|) . Since both D"^ and Fi 
can be expressed purely in terms of scalar fields they can be used to write the complete 
scalar potential as 

F(0*, </)) = f:f^ + ^y1 = -w:w' + Iy1 9I • (6-12) 

a a 
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1 






Table 6.1: The interaction states of MSSM and the respective gauge charges. There are 
also three generations of (s) quarks and (s)leptons. 

6.1.2 Minimal Supersymmetric Standard Model 

The Minimal Supersymmetric Standard Model (MSSM) is the model that contains all 
of the possible = 1 SUSY interactions and special cases of this model are generally 
analyzed phenomenologically. This model has one superpartner for every SM particle 
and an extra SU{2)l Higgs doublet. The extra Higgs doublet is due to the requirement 
that the superpotential be an analytic function. This prevents the addition of a Yukawa 
term like Q(j)^UR and instead a second Higgs doublet is needed. Table 16.11 shows the 
particle content in terms of the interaction eigenstates. We also introduce the convention 
that all superpartners of fermions are names with the letter s preceding the name and 
all superpartners of bosons have ino appended to the name. For example the partner 
of an electron is a selectron and the partner of a B is a Bino. Also, by convention the 
symbol for a superpartner field is to put a tilde on top of the field symbol, e.g. 5^ B. 
Bosons with a superscript /i are vector bosons, whereas those without are scalar bosons. 

Just like in the SM, several of the interaction states mix to form the mass states. 
Table shows which interaction states mix to form which mass states; the mass states 
are denoted by the word for the physical particles that are expected to be observed. The 
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Mixing Fermions 



Mass States 



Mixing Bosons 



Mass States 



quarks 
electrons, neutrinos 
neutralinos, x° 
charginos, xf 



Ql, ur, da 



Hi, H2 
9^ 



squarks 
selectrons,sneutrino 
Z^, photon 

CP-even/odd Higgs, Charged Higgs 
gluon 



9 



gluino 



Table 6.2: The interaction states that mix and yield the mass states. The mass states 
are given by name. 

mixture of the (s) quarks given in the table is to signify that the left and right handed 
(s)quarks mix, not the up- and down-type (s)quarks. The Higgs sector isn't quite as 
straightforward. Each Higgs doublet contains a charged part and a neutral part. The 
real parts of the two neutral Higgs doublets mix to form the CP even Higgs states, often 
denoted by h and H. Here the lower case is the lighter Higgs, the one that is expected 
in the SM. The imaginary parts of the neutral Higgs doublets mix to form the CP odd 
Higgs, A. The remaining charged parts of both doublets then mix to form two charged 
Higgs bosons. Also, often the convention if„ and is used [90]. By the definitions 
given here this is equivalent to Hi = H^ and H2 = H^^. The mass states of the four 
neutralinos are denoted by Xi,2,3,4 order of mass. The four mass state charginos are 
denoted by Xi,2- 

6.1.2.1 Soft Breaking 

As none of the superpartners of the SM particles have been observed this implies that 
supersymmetry must be a broken symmetry. If we refer back to the fine tuning argument 
we see that in order for the quadratically divergent parts to still cancel, the dimensionless 
couplings must cancel (i.e. A5 = I'^/l^)- This leads us to only consider "soft" breaking 
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^ — -^SUSY + -^soft (6.13) 

where >Csusy are all the terms preserving supersymmetry and >Csoft a-re all the terms that 
softly break supersymmetry. To ensure renormalizability and to maintain the natural 
cancellation of quadratic divergences, >Csoft must have only mass terms and couplings 
with positive mass dimension. 

It is the existence of the soft term that allows all of the superpartners of the SM 
particles to be heavier than the top quark. In the absence of electroweak symmetry 
breaking, all of the SM fields would be massless. This isn't true for the superpartners. 
The scalars can have a mass term in the Lagrangian of the form The gauginos 

and Higgsinos also do not require electroweak breaking in order to acquire a mass due 
to the fact that they are fermions in a real representation of their gauge group. 

If the largest mass term occurring in jCgoft is TTigoit- We can see that the corrections 
to the Higgs mass-squared, due to the SUSY particles, must vanish as mgoft 0. This 
means that these corrections cannot be quadratic in Auv and the corrections must be 
of the form 

^ (6.14) 



In (Ac/y/msoft) + • • • 



167r2 

If we take A of the order 1 and A^y of order Mp, we find that the lightest SUSY particles 
should be about 1 TeV, as stated previously. This is what leads to the optimism that 
SUSY will be discovered at the LHC. 

The actual soft breaking terms added to MSSM are 

£,oft = (M;,At'^A" + c.c.) - {m^y. 0*^0, - 06^^0,0,- + ^a^'^^k + c.c}j . (6.15) 
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There is another term that could be added without losing the renormalizability of the 
theory. This term is — |c^^0**0j0fc + c.c. but this is often not included as it can lead to 
quadratic divergences in some models. It is also evident that the terms in do break 

supersymmetry because they are not expressions of the supermultiplets. This is because 
mass terms for the fermions can be reabsorbed into a redefinition of the superpotential 
and the (m^)* and terms. 

To complete the definition of the MSSM we must also specify the superpotential. 
This is given as 

W^MSSM = Y^'^uH.Q - YfflH^Q - Y^'jcH^i + fxH.H^. (6.16) 

Here the fields Q, i, Hi, H2, u, d, e are the chiral supermultiplets, not just the SM fields. 
We can see the family Yukawa matrices and the Higgs fields which are used to give the 
masses. We also see the supersymmetric equivalent to the /i term in the SM. 

This section has been a simple introduction to SUSY models, SUSY Lagrangians and 
MSSM. For a more detailed discussion of the development of all these topics see [90-92]. 
Next the concept of the effective potential and its uses is given and finally the software 
Effective is discussed which uses the effective potential to study = 1 SUSY models. 

6.2 Effective Potential 

The effective potential was originally introduced by Euler and Heisenberg [97] and further 
expanded by Schwinger [98] . This was later applied to studies of spontaneous symmetry 
breaking by Goldstone, Salam and Weinberg [99]. The development of the effective 
potential given in this section, as well as a complete summary of the effective potential 
and its uses is given in [100]. 
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The effective potential is capable of providing quite a lot of information about a 
theory, while working with a simpler expression. Here is presented the development of 
the effective potential and in subsequent sections the various things that can be derived 
from an effective potential are given. These are all capabihties of the software, Effective, 
which is discussed in the last section. 

We start by giving the theoretical derivation of the effective potential. With this 
derivation it can be shown that the the one-loop contribution to the effective potential 
can be derived in a model independent way. Unfortunately, the higher order corrections 
cannot be developed in a model independent way. Because of the model dependence, 
the two-loop effective potential is not implemented in Effective, though future extensions 
of the code could provide this functionality. 

6.2.1 Generating Functionals 

As an example we start with a theory described by the scalar field with a Lagrangian 
density jC{(p{x)}. The action is then 




(6.17) 



The generating functional is the vacuum-to- vacuum expectation value {Oout\ Oin)j and 
is given by the path-integral representation. 




(6.18) 



where 




(6.19) 
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Using ()6.18|1 we define 

Z[j]=exp{zW[j]}, (6.20) 

where W [j] is known as the connected generating functional. The effective action T [0] 
is just the Legendre transformation of ()6.20|) and is given by 

r [0] = W[j] - I d^'x^-^jix), (6.21) 



where 



6W \j] 

m = -TTT- (6.22) 

6j{x) 



4>{x) is the weighted average of the fiuctuations of the field. In a translationally 
invariant theory, which are the ones Effective is designed to deal with, (j){x) is a constant 



0(x) = (6.23) 
The effective potential can then be defined as 




(6.24) 



which can be written as an expansion as 

oo _ 

Kff (0c) = -Yl -f'^^r^") {P^ = 0) , (6.25) 

n=0 

where F^") are the one-particle irreducible (IPI) Green functions. Minimizing the effec- 
tive potential over the constant fields, (pc gives the vacuum state of the theory [3]. 
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The tree-level effective potential is identical to the classical effective potential. This is 
simply 



V, 



tree 



c 



(6.26) 



The one-loop contribution can be written in closed form for any theory containing fields 
of spin 0, ^, or 1. 

We show here the one-loop correction for a model with one self-interacting scalar 
field described by the Lagrangian 



C = -d^(j)d,(j)-Vo{(j)), 



(6.27) 



where Vq is the tree-level potential given by 



(6.28) 



As we expressed in ()6.25j) . the one- loop correction to the tree-level effective potential 
is given by the sum of all IPI diagrams with a single loop and zero external momenta. 
The n-th diagram has n propagators, n vertices and 2n external legs. The n propagators 
contribute a factor of i^{p'^ — iv? + it)^'^ , as we saw in chapter ^ Each pair of the external 
lines contributes a factor of 0^" and each vertex a factor of —i A/2. There is also a global 
symmetry factor of l/2n. 



This gives 



d^p 1 

^ (27r)4 2^ 

n=l " ^ ' 

i f d'^p 
"2 J 



— + ie 

p2 _ ^2 _|_ 



(6.29) 
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In fl(i.2()|l we did not include a normalization factor of exp{iiy[0]}. This is needed to 
correctly generate ()(j.29j) . Without it a shift is introduced to the effective potential. 
After a Wick rotation and using the unnormalized generating functional we find 



log [p^ + m^(0c)] 



(6.30) 



where the momentum is now in Euclidean space and the mass is the shifted mass given 
by 



^ ('^c) = T75 — • (6.31) 



This is finally 



1 4..^^ ^'i<l^c) 3 

-m (0c) In 



(6.32) 



in the DR renormalization scheme [101] where fi is the renormalization scale. This result 
can be trivially generalized to the case of A^^ complex scalar fields each described by the 
Lagrangian, 

c = d^rd,<i>i-vo{rAi)- (6.33) 

The one-loop contribution in the DR renormalization scheme is 



647r2 



Tr 



m: 



In 



m: 



2f J^a 



(6.34) 



where Mg is the mass matrix of the fields given by 



d^V 
d(i)id4>^ 



(6.35) 



Eq. ()6.34p can be generalized to fermions obeying the Dirac equation and gauge 
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bosons [100, 102] as 



E^(2.+ l)(-ir^'TV.[M^(,".,l)(.„M)_3) 



(6.36) 



where the sum is over the spins, 0, 2, and 1, and Tr^ is the trace is over the mass matrix 



of spin, s. Defining the "Supertrace" as 



STr/(X) = ^(25, + 1)(-1)2^7(X,), 



(6.37) 



which is the spin-weighted trace, we have the one-loop contribution in the DR renor- 
mahzation scheme as 



The fields defined as representations of the groups in a field theory are interaction eigen- 
states of the theory. These are states in which the interactions are directly governed by 
the respective couphngs. The particles that are observed are the mass eigenstates of the 
theory. These states are a mixture of interaction eigenstates, and thus the interactions 
of these states are a mixture of the couplings in the theory. The interaction states that 
mix to form the mass states, as well as the mixing angles, can be derived from the mass 
matrices. In should be noted that the mass matrices are not derivable from the effective 
potential. Rather, the value of the VEVs chosen from the minimization of the effec- 
tive potential determine the mass spectrum of the particles. Derivation of the masses 
requires the full potential. 




(6.38) 



6.3 Mass Matrices 
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For complex scalars the mass matrix is defined as 

d^V 



2\« 



(6.39) 



4)i{x)^4>i, 



where is the expectation value of the the scalar field (f)i{x). Similarly for vectors the 
mass matrix is 



(6.40) 



4>i{x)^4ii, 



dAadA^ 

The mass matrix for fermions is not defined as a mass squared matrix. Instead for 
complex fermions it is given by 



d^V 



(6.41) 

4>i{x)^4iic 



Diagonalizing these mass matrices gives 

Mly = U^DlyV, Mj = U^DjV, (6.42) 

where U and V are unitary matrices which define the mixing angles of the interaction 
states. D is the diagonal matrix whose elements are the masses of the mass eigenstates. 

Eqs. (16.39116.411) don't indicate whether the potential is to be taken at tree-level, one- 
loop or some other order. That is because the masses also have higher order corrections. 
If the tree-level potential is used, it yields the tree-level mass matrices. Likewise, if the 
one-loop potential is used it yields the one-loop mass matrix. For the exact mass term 
the exact potential would be needed. 

As we saw earlier, the one-loop effective potential can be written in terms of the mass 
matrices only. The matrices in ()6.38|1 are the tree level matrices. From 1)6.39116.4111 we 
can see that the one-loop masses are derived from the full one-loop potential. This would 
require the calculation of the full one-loop potential, which is much more complicated. 
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Instead the one-loop corrections to the mass matrices can be computed in another way. 
If we refer back to chapter the masses of a field are related to the two-point Green's 
function of that field. In fact, it can be shown that the one- loop correction to the mass 
is just the one-loop correction to the two-point Green's function [3]. 

Before we present the one-loop corrections to the mass matrices, the idea of a ghost 
field must be explained. We saw that in chapter ^ that the propagator for a gauge 
boson is given as ^j^. This is not quite true for non-Abelian groups. Instead using the 
Faddeev-Popov method to quantize the field and the gauge condition 



G{A) = d^Alix) - u^ix) = 0, 



(6.43) 



where cu" is a Gaussian weight, we find the propagator is 



[A',ix)A:ix)) 



+ ie 



(6.44) 



The Faddeev-Popov quantization inserts a determinant of the form det y — I where 



( 5G{A°'y 

a is now used to represent the infinitesimal form of the gauge transformation. This 
determinant is zero for an Abelian group but non-zero for a non-Abelian group. Instead 
Faddeev and Popov showed it can be expressed as 



, f6GiA 



det ( -d^'D 
9 



VcVc exp 



t I d^xc{-d^'D^)c 



(6.45) 



where c and c are now anti-commuting fields which are scalars under Lorentz trans- 
formations. These new fields are known as ghost fields and introduce new Feynman 
rules that must be considered to ensure that the value of an S'-matrix calculation is 
gauge-invariant . 
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For a spontaneously broken symmetry, we impose the gauge condition 

G{A) = -^{d.A^ - i9v<i>) = 0, (6.46) 

where is the Goldstone boson associated with the breaking and v is the VEV of the 
field which leads to the breaking. Using this form of the gauge fixing condition gives the 
massive gauge boson the propagator 



(A^Ai,) 



p2 _ ^2 



9t,u - — — 1 -C) 



(6.47) 



which is the same value for the gauge boson propagator shown before for the massless 
case. This also gives the Goldstone boson, (p, a mass of ^m^ where uia is the mass of 
the gauge boson. The ghost also acquires the same mass, ^m^, and this must be taken 
into account when computing Feynman rules of a process. 

The one-loop corrections for a field can be given in a gauge independent form by 
explicitly keeping ^. That has been done for the one- loop corrections given here. In 
Effective, the ghost contributions are hard to automate. Instead we work in the Landau 
gauge, ^ ^ where the masses of the Goldstone bosons remain zero and the ghost terms 
do not contribute to the calculations. 

The one-loop corrections can be written in terms of the Passarino-Veltman [103-105] 
functions. These are defined by 



Ao{m'';f/) = 1677^^ 
So(p^m?,m^;//^) = 167rV"'^ 



id'^q 



idf-q 1 



(27r)'^ \q^ — ml + ie] [{q — p)^ — ml + ie] ' 
id'^q q^ 



{2tiY [9^ — m\ + ie] [{q — — ml + ie] ' 
P,P.B,,{p',mlml-^i') + g,^Boo{p\ml,ml-^i') (6.48) 
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id'^q q q„ 



{2nY [g2 — rnl + ie] [{q — pY — m| + ie] ' 



where is the renormahzation scale and the integrals are regularized in d = 4 — 2e 
dimensions. The first two of these can be integrated and expressed (expanding to 0{e)) 



as 



A^{m'■^i') = m^Q-l + In^^+C»(e), (6.49) 
So(p^m2,m2;//2) = ^ _ in ^ - /5(x+) - /b(x_)z + 0(e), (6.50) 

where l/e = 1/e — 7B + ln 47r, 

x± = , /b(x) = ln(l -x) -xln (1 - - 1 - 1, (6.51) 

and s — p'^ — ml + ml. The other Passarino-Veltman functions can be decomposed 
into the two scalar functions, Aq and Bq [104]. There are more functions for higher 
loop corrections that are not given here. These can be decomposed into their respective 
higher-loop scalar functions (e.g. Co, -Do, ■ ■ ■ )■ 

The corrections are summed over all internal fields and the couplings are given from 
the Lagrangian. There are two ways that the one-loop corrections can be applied. 
First one could find the diagonal form of the tree-level matrix. This defines the mass 
eigenstates and therefore the couphngs of these states. The corrections can then be 
directly applied to the tree-level masses given the couplings defined in this manner. This 
is not exactly the one-loop correction, however. This will give a close approximation to 
the one-loop masses, but will only give the tree-level mixings. If instead the one-loop 
corrections are applied to the undiagonalized matrix, with corrections to the off-diagonal 
elements, then the diagonalization of this corrected matrix will give the correct one-loop 



■^Not to be confused with jj, for the Lorentz index 
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/2 V2 5 

(d) (e) (f) 



Figure 6.2: The different diagrams that contribute to the one-loop correction of the 
scalar mass matrix element {M^)ab- 

masses and mixings. Both techniques have been implemented in Effective, but the latter 
is much more time intensive and if only the masses are desired, applying the corrections 
to the mass eigenstates is a good approximation. 



6.3.1 Scalar One-Loop Corrections 

We first give the one-loop corrections to the scalar masses. Figure 16.21 shows the cor- 
rections to the scalar mass matrix term (M^)|J. These diagrams are summed over all 
possible intermediate states given in the theory. The results given here have been cal- 
culated with the help of the computer software FeynCalc [106] and FormCalc [107]. 

The form for the contribution of the scalar four-point vertex, given by the first 
diagram in figure 16.21 is 

{^^-%{v^m';i^') = %-^Ao{m^;^.'), (6.52) 
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where is the four-point couphng of the fields a and b with an arbitrary scalar field 

s. The second diagram in figure l(j.2l is the contribution of the four-point vertex given 
with a vector loop. This correction is 

^j^iavvb)y^ m'; = ^^^^ {2m' - 3Ao{m') + ^Ao{^m')) , (6.53) 

The diagram given by fig. I(i.2b is of the scalar three-point vertex. This is given by 

(S(— ^));(p^m?,m^;/.^) = -^—^ i?o(p^m?,m^) (6.54) 

where C^^^'^^^^ is the coupling of field a to the two fields in the loop and C^^^^'-^'''^ is the 
coupling of field b. 

There is no fermion four-point diagram as this type of term would have a mass 
dimension greater than four in the Lagrangian. Instead there is only a three-point 
correction. This term (fig. I6.3l i) is given by the expression 

{j:(^f^f^%{p',ml,ml-f.') = {Ao(m?;/i^) + Ao(m^;/.2) + 

((mi +m2)^ - p'')Bo{p'^,ml,ml; ij!^)} . (6.55) 

The vector three-point term, given by the fifth diagram in figure 16. 2| is 

ry{aViV2)r<(yiy2b) r r? r? 

- «) si - i^Bi' + 

7712 ^2 V "^2 / ^2 

2 

(^00 - ^00 + ^00 - ^oo) + 



2 2 

mfm2 



{bIi - + - B\l) + 



mfm2 



2e^(i?i-5P)V (6.56) 

I7l2 ) 
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where 



Al = Ao{ml;f?), Al = Ao{ml; fi^) , A^ = Ao{^ml; fi^) , A^^ = Aoi^mj; fx^) , 
Bl = B,^{p^, ml, ml; /i^), El = B^{p^, ml, ^mf; /x^), = B^{p^, ml, ^ml; /x^), 
B'^ ^ B,ip',^ml,^ml;ix'). (6.57) 

Remarkably, in the 't Hooft-Feynman gauge (.^ ^ 1) this reduces to 

(E^'^v.v.^ip^mlml^fx') = ^^^^ (l - 2i?o(p^ m?, m^; /x^)) . (6.58) 

The last one- loop corrections to the scalar mass matrix, (M^)^, is given by the last 
diagram in figure 16.21 This is the correction given by the vector-scalar three-point 
vertex, given by 

(E("^^^)):(p^m^,m^;;.2) = {{p' - m^A^ - {p^ - - ml)Al+ 

{mUp' + ml)~{p'-mlf) 5° - 
(p2 - ml^B] - 2m^p2^°} , (6.59) 

where the vector corresponds to loop field 1 in the previous abbreviations and the scalar 
field is loop field 2. 



6.3.2 Vector One-Loop Corrections 



The one-loop corrections to the vector fields can also be generated using the same two 
programs FeynCalc and FormCalc. Figure shows the diagrams which give corrections 
to the vector mass. 
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Sl 
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a \ 



S2 




(c) 




(d) 




(e) 



(f) 



G 



G 
(g) 



Figure 6.3: The different diagrams that contribute to the one-loop correction of the 
vector mass matrix element {My)ab- 

The first diagram in figure 16.31 leads to a correction of the form 



( 



KT'^'X {P\<rnl;f,') = ^ B,,{p\ ml ml [^g,^ - ^ ) • (6.60) 



The contribution from the fermion loop is gauge independent and is given by 



{^t^'^^X{vlmlml^^' 



(7(a/l/2)(7{/l/2fe) 



{Ao{ml fi^) + Ao{ml + 



((mi - ms)^ - p^) Bq{p^, ml ml /i^ 



45oo(/,mi,m2;/i^)} ( g^,,, 



p2 



(6.61) 



and fig. 16.3b shows the scalar four-point correction. This is 



iKV'Xip^mlrnl^, 



-4o(ms;/i)U^^ 



327r2 



pi 



(6.62) 
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The next two diagrams are of the four-point vector loop and the three-point vector- 
scalar loop. These are given by 

Kr'X (P^ /^^) = ^^^^ (2m^ - 5Ao(m2; /i^) - ^Ao{^m'; /i^)) [g,. - ^ j , 

(6.63) 

and 

„-/^(ays)(^{Vsfe) / n n \ 

fTj(aVsb)\'^ / 2 _2 2. 2\ _ '''^ (2 r2 n2 , r12\ / ^ /^M/^^ \ 

li^A.:. ')f,{P ,mv,^s,l^ ) - 32^2^2 l"^V^O - ^00 + ^00 j " I , 

(6.64) 

where we have used the same abbreviations ()6.57|1 and the vector corresponds to field 1 
in the notation and the scalar is field 2. The remaining two diagrams are for the vector 
three-point correction and the ghost terms. The vector three-point correction is quite 
complicated and is given by 

r("-viV2)r<(yiy2b) co 

(P2 - + 0^0 - (P2 - - 0^0 + 

{ml{l - Ki) + p2(2/«i + 5 - p2)5° + 
{ml{K, + e)+p\p2-2K,))Bl + 
{ki + 11 - 2p2 - 5pi + piP2)S°o + 

(2p2 - Ki - PiP2)BIq + (5pi - pip2 - 0^00 + P1P2BII + 
-2p^Bl + Aip'Bl + 

4pi/ {Bl, - Bl) } [g,^ - ^) , (6.65) 
where the abbreviations (16.571) are used and = pi = ^ and P2 = In the 

^^^^^^^^^^ TTT.^ tFL-^ tFht^ 

't Hooft-Feynman gauge this reduces to 

r'(aViy2)r'(^i^2fe) 

{Ii^^^^y^X(p\n^ln^l^^-) = ^-^^ - 2(m? + m^) + 
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^ V 




(a) (b) 

Figure 6.4: The different diagrams that contribute to the one-loop correction of the 
fermion mass matrix element (M/)^. 

2 {mj + 3/) 5o(p^ mj, ml; - 2KiAo{ml; + 
(10 + pl)5oo(p^ ml, ml; /i^) + 

2p'B, (p^ ml ml /i^) } [g,^ - ^) . (6.66) 

It must also be noted that the gauge bosons from an unbroken group are massless. These 
bosons only couple to themselves and the correction is proportional to p^. The correction 
is to be evaluated at = m^, which for bosons from an unbroken group is zero. 

The remaining one-loop correction to the vector mass matrix is from the ghost terms. 
This is the last diagram in figure 16.31 and is given by 

(n^r^O: (^'' -c; /^^) = ^^^^^^^Booipl ^ml, iml; /.^) [g,. - ^) , (6.67) 
from which it is easy to see that the ghost contribution is zero in the Landau gauge. 



6.3.3 Fermion One-Loop Corrections 

There are only two mass corrections to the fermion masses. The diagrams of these are 
given in fig. 16.41 The first one is the correction due to a scalar-fermion three-point 
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coupling. This is 

(S(-/'')); (p^ ml mj; /i^) = ^^^^ ^B,{p\ ml mj; (6.68) 

The mass correction is given by bm = S(^)|^=m- "^^^^ approximated by 5^ ~ 



Lastly, the correction given by fig. 16.4b is from a vector boson three-point coupling. 
This is 

(S(«^/^)): {plml,m};^^^) = ^-^^ ^ [fi? - (1 - 0(^0^0 + , (6-69) 

where the abbreviations given earlier are used and the index 1 refers to the vector and 
the index 2 is the fermion. 



6.4 Renormalization Group Equations 

The basic concept of renormalization stems from the observation that the divergences 
in the one-loop graphs amount to shifts in the parameters of the action. For example, 
as we saw in the last section they change the mass of a field. Renormalization is the 
procedure of cancelling the divergences from these shifts by introducing counterterms 
into the Lagrangian. These counter terms are defined such that they exactly cancel the 
divergent quantity in the one- loop correction but leave the finite parts untouched [5]. 

We can first consider the 0^ theory. The Lagrangian is given by 

C = (90o) V2 - ml4/2 - (7o0o/4! (6.70) 



We can then rescale the field by (pQ = Z^^'^cp, so that in terms of the new 'renormalized 
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field' (f) we have 

C = Z{d(f))y2 - Zml(t?/2 - g^Z"^ (j)'^ / (6.71) 
We can then write this in terms of the counterterms 

C = (90)V2 - m20V2 - 

+5z{d(t)f/2 - 5m^(t?/2 - dgct)^^] (6.72) 

We can see that 5z = Z — 1, = Zm^ — and 5g = Z'^g^ — g. These counterterms 
are adjusted to cancel the divergences of each term. 

To renormalize a theory renormalization conditions must be imposed. The general 
conditions that are imposed are to set the contributions from the IPI diagrams to the 
two-point Green's functions and the derivatives of these diagrams to zero at a spacelike 
momentum, = — /i^. /i is then known as the renormalization scale. The conditions are 
imposed at a spacelike momentum in order to avoid threshold singularities. Similarly, 
the IPI contributions to the three- and four-point Green's functions are defined such 
that the coupling at that scale is exactly the physical coupling (e.g. g from ()6.72|) ). 

The renormalization scale is completely arbitrary. The theory can just as easily be 
defined at another renormalization scale. This implies that the theory and the redefined 
couplings should depend on the scale in such a way that physical calculations are in- 
dependent of the scale. For example if we take an ra-point Green's function we require 

^Q{n) ^ ^^^Sfi + ^-^5A = n6r]G^''\ (6.73) 
Ofi oX 

This is for either a massless field or if the mass of the field is considered as on of the 
vertices that needs to be renormalized. 6X is the counterterm corresponding to the 
coupling A and 6r] is the rescaling of the external field in the n-point Green's function. 
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Defining 



^ dX drj 



we can now give the Callan-Symanzik equation [108, 109]. This is 



(6.74) 



d ^ d 



G^""^ = 0. 



(6.75) 



The P and 7 functions are the same between different Green's functions. Knowing 
this we can find the 7 function from the 2-point Green's function of the field. The 2-point 
Green's function can be expressed in terms of the loop diagrams and the counterterms. 
Doing so 7 can be shown to be 

1 f) 

(6.76) 



1 



The counterterm can be written in the form (at lowest order) 



A2 

5z = ^ log — + finite, 



(6.77) 



in order to cancel the divergent logarithm in G^^^ . This means that we can find 7 simply 
by determining A, the coefficient of the divergent piece of the IPI diagrams. Since 7 is 
related to the field renormalization it is known as the anomalous dimension of the field. 

Knowing the anomalous dimensions, the /3 functions can be derived by using the 
appropriate Green's function (G^"^ for an n-point coupling). The n-point Green's func- 
tion is given by the sum of the tree-level diagram, the IPI loop diagrams of the vertex, 
the vertex counterterm and the external leg corrections. This can be expressed for the 
n-point coupling at one- loop order, as 



-ig - iB log — ^ - i5g + {-ig) ^ ( A^ log — ^ - 5^, J 

P 4 \ Pi / 



, (6.78) 
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plus the finite terms. Applying the Callan-Symanzik equation we find 

d 



We can see that the second term on the right is just the sum of the anomalous dimensions 
of all the external particles in the Green's function. At one-loop, 6g is given by an 
expression of the form 

A2 

5g = -B log — + finite. (6.80) 

This means that we just have to compute the divergent pieces of the IPI corrections to 
the coupling and Pg is given by 



f3g = -2B + gY,l^■ (6.81) 



Most theories have many couplings and fields associated with them. This requires 
the procedure that produced ()6.75|) to be applied to a set of couplings, A/, and a set of 
fields 0i. 

/i^ + /5/^ + 7i0i^ O{Xi,mi-fi)=0, (6.82) 

We saw earlier that using the two-, three- and four-point Green's functions, eq. fl6.82|l can 
be applied to determine the form of all jSj and all 7^, known as the renormalization group 
equations [110-112]. This technique requires the calculation of several IPI diagrams 
and is quite computationally intensive. Instead there are some couplings which can be 
computed either completely for an arbitrary theory, or by applying the Callan-Symanzik 
equation to a different set of observables. 

The running of the gauge group couplings can be computed from many different 
vertices. They all require the boson self-energy which is given by the diagrams in fig. 
16.31 The divergences of these diagrams define the anomalous dimension of the gauge 
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boson field. These sum to give 



7b = 



(47r)2 




(6.83) 



where i is the sum over fermions, j is the sum over scalars, C2{G) = f°-^'^f°-'^'^ and 
C{r) = Tr[t"t^]- O'^^ can then use the fermion-fermion- vector coupling, as in [3], or one 
could just as easily use the gauge boson three- or four-point couplings (in non-Abelian 
theories). The final result is gauge independent and given in general by the expression 



where again the index i indicates the fermions and j the scalars. For the case of Weyl 



The SM has Yukawa couplings, gauge couplings and couplings of the Higgs boson 
(A and /i). These are all dependent on the scale ^, Q, at which they are evaluated. A 
common argument for beyond standard model (BSM) physics is that given the particles 
in the SM, the couplings, a, a^eak and Og, will never meet at a single ultimate scale, 
known as the unification scale. Though there is no proof that indicates these should 
unify, it is believed to be a desired and expected result of a theory. Figure 16.51 shows 
the running of the different a's for the particles in the EW theory with one generation 
(black, sohd), the SM with three generations (red, dashed) and the MSSM with three 
generations (blue, dot-dashed). We can see from this figure that the couplings from the 
MSSM all meet at a single scale. The derivation of the f3 functions and the running 
of the couplings in this figure were all done by entering the models into Effective. In 
unification theories, the U{1) coupling is multiplied by ■\/5/3 and this is also done here. 

'^Here Q is used instead of /i to avoid confusion with the Higgs couphng, fj,. 




(6.84) 



spinors, the sum of the fermions has an additional factor of a half (e.g. 
|), and similarly for real scalars. 



I rather than 
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Running couplings 




10^ 10'' 10'^ 10« 10'^ 10^^ lO^'' 10^^ 10^« 

Q/GeV 

Figure 6.5: The scale dependence of the different gauge couphngs. The black (solid) 
lines are for the EW model with one fermion generation. The red (dashed) lines are for 
the SM with all three fermion generations. The blue (dot-dashed) lines show the the 
couplings for the set of particles in the MSSM. 

In Effective the observable that is used to define some of the couplings is the effective 
potential [100,113]. The one-loop effective potential Vo + Vi, where Vi is given in 



can be used to define the RGEs for all the couplings that appear in the potential. This 
is the set of couplings for the fields that may acquire a non-zero VEV. If we apply the 
Callan-Symanzik equation to the one-loop effective potential, where now the derivatives 
with respect to the fields are replaced by derivatives with respect to the VEVs, this 
yields 



d 



V. 



327r2 



STr 



(6.85) 



Using this one can compare the powers of the VEVs (0jc) on the left hand side of the 
equation to those of the right and derive the relations of the (3 functions, given the 
anomalous dimensions, 7i. 

To compute the remaining (3 functions, we must use (16.811) . This requires us to 
compute the IPI diagrams for the three- and four-point vertices given in fig. 16.61 Since 
we have already derived the {3 functions of the gauge group couplings we do not need to 
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Figure 6.6: The three- and four-point diagrams that must have their IPI diagrams 
computed to define the complete set of /? functions for a theory. 

consider the vertices with a gauge boson as an external leg. 



6.5 Effective 



Effective is a computer program that is able to generate the effective potential, the mass 
spectrum and the RGEs of a theory. It is originally based on the work in [114] and uses 
the GiNaC [115] computer algebra system to do algebraic manipulations. This package 
provides for evaluations and simplifications of complex algebraic structures within the 
C-I-+ environment. As this is a C-I--I- library, it also allows for expansion. New functions 
and objects can be created that can then have their own set of functions that are 
used during simplifications and evaluations. It is the ability to integrate the computer 
algebra with the standard C-I--I- functionality and to be able to define new objects that 
interact with the computer algebra package that make GiNaC an ideal library for building 
Effective on. 
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My contribution to Effective covers a broad spectrum. I have restructured and written 
the code to work with the algebra package, GiNaC. The original one-loop mass corrections 
in the code were only for scalar particles. These were taken from the literature and did 
not all use the same gauge and as a complete correction, were wrong. As has been 
seen earlier, these have been rederived in an arbitrary gauge for all different spins and 
implemented in Effective in the Landau gauge. Also, I have done the development and 
implementation of the RGEs in Effective. 

Effective has currently been designed to work with N — 1 SUSY. Effective itself is a 
package that can be used to define a SUSY model. Once this model is defined a user 
can interact with it. It is these different interactions that are discussed in this section. 

6.5.1 Model Definition 

A model in Effective is defined by a set of gauge groups, fields and additional interactions. 
These are then used to define a Lagrangian and an effective potential. The Lagrangian 
can be used to define the couplings of the theory and the mass matrices. As wc saw in 
the last section, the effective potential as well as the Green's functions can define the 
RGEs of the theory. 

A gauge group is just that. In Effective this is a class that defines the structure 
constants and the generators. This also defines the dimension of the indices in the 
fundamental and adjoint representations. It is also possible to have several gauge groups 
of the same group structure defined. For example if a user wants to test a theory with 
two SU{2) groups this can be done. Each gauge group has an extra flag to indicate 
the 'line' that it is part of. If someone wants to use two gauges with the same group 
structure, they just have to deflne a new line for the new group. When indices are being 
contracted, these lines are compared and only when they are identical will a contraction 
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occur. Effective currently has the SU{3), SU{2) and U{1) groups implemented, but it is 
easy to expand and add new groups. An example of this is given in the appendix. 

Once the groups have been defined, the gauge fields can be added. These are the 
fields which mediate the interactions. It is possible to create gauge bosons and their 
superpartners, gauge fermions. When the fields are created the relevant terms can be 
added to the Lagrangian. This is given in section EH] but is repeated here for convenience 

/:gauge = --^F'^^'^F;, - ^AtV^D,,A^ (6.86) 

where A'^ is the gauge fermion and F^^ is the standard tensor formed by a gauge boson. 
In = 1 SUSY there are no other terms which depend only on the gauge mediators. It 
is possible to extend Effective to provide for 7^ 1 SUSY. If this were done then a new 
set of terms will also need to be added to the Lagrangian. 

After the gauge groups and the gauge fields are defined the matter fields can be 
given. These are fields that can have charges under more than one group. These fields 
also define the other gauge invariant terms. For a given fermion, ip, and superpartner 
scalar, 0, we have the terms 

/:mattcr = (T^ D + (^^0)^ - V2g, [(0*T» A"^ + C.C.] + ]^gf {<j)*T''<f)f , (6.87) 

where A is again the gauge mediating fermion, is the covariant derivative for the 
groups the fields are charged under, and gi is the coupling of the i-th gauge group. This 
expression is summed over all the gauge groups of the matter fields. It is also at this 
point that the VEVs of the scalar fields are defined. A set of parameters can be defined 
and combinations of these parameters can be given as the VEVs. For example, in MSSM 
one often works with the VEVs in terms of v and tan (3, not the actual VEVs, Vi and 
V2- This can also be done in Effective so the set of parameters that one works with 



CHAPTER 6. EFFECTIVE POTENTIAL ANALYSIS: EFFECTIVE 205 

is consistent with the set of parameters given in the hterature. An example of this is 
given in the appendix. Unfortunately, some of the routines in Effective require the actual 
VEVs to be entered. These include the RGEs and the minimization of the potential. If 
neither of these functions is needed, however, then the set of parameters used in a paper 
can also be used in Effective. 

Once the fields of the model are given, the superpotential can be defined. In Effective 
this is not given as a function of superfields, as the definition of the superpotential is, but 
the user must specify this in terms of the relevant fields (scalars, fermions, etc.). This 
is a feature of the code that could be improved to provide the appropriate superfield 
formalism. Once the superpotential is defined then the following terms are added to the 
Lagrangian 

jC^ = {W'^jjjijjjj + c.c.) - W*W\ (6.88) 

where we saw before that 



d^W dW 



After the definition of the superpotential the last remaining ingredient of a model 
is the symmetry breaking terms. These are usually the soft breaking terms of a SUSY 
model but Effective can be used to work with non-SUSY models as well. In this type of 
model one would define any other gauge invariant terms that are desired for the model 
but have not been added up to this point (e.g. the Yukawa terms). 



6.5.2 Mass Matrices 



Once the VEVs have been defined the code can generate the tree-level mass matrices. 
Effective is able to do this and produce the matrices in analytic form. It is important to 



CHAPTER 6. EFFECTIVE POTENTIAL ANALYSIS: EFFECTIVE 206 

note that in Effective the real and imaginary parts of the fields are treated independently. 
The separation of these parts has been implemented to provide a mechanism for studying 
CP violating models. This separation causes several complications when comparing to 
the literature though. 

Often in the literature the states are given as the mass eigenstates. These are a 
combination of real and imaginary parts of the the interaction states. It is these real 
and imaginary combinations that must be compared to the results of Effective for a true 
comparison. Often this is a tedious task. Instead it is usually more convincing, and 
easier, to check the numerical results of Effective against the literature. 

Another complication that arises when implementing formulae in Effective is that 
Weyl spinors are used, not Dirac ones which are often used in literature. This affects the 
degrees of freedom of the fields in Effective. Usually one considers a fermionic field to 
have a spin factor of —2 from the spin coefficient (2s + 1)(— 1)^'^. Instead since Effective 
uses Weyl spinors, this has a factor of | attached to it, giving —1. The difference between 
Dirac spinors and Weyl spinors must also be taken into account when implementing the 
one-loop corrections given in section 

As mentioned earlier. Effective uses the Landau gauge = 0) to avoid the need for 
ghost contributions. The one-loop corrections to the mass matrices have been imple- 
mented in Effective in this gauge. These corrections are given in numerical form only. 
There are two ways in which the corrections can be applied. The first is that the tree- 
level matrices are diagonalized and the mixing matrices are set. The corrections are then 
applied directly to the mass eigenstates. This is only an approximation to the one-loop 
corrections and provides only the tree-level mixing matrices. The other way in which the 
one-loop corrections can be applied is to add the one-loop corrections of the interaction 
eigenstates to the tree-level mass matrix. This new matrix is then diagonalized and this 
process defines the one-loop mixing matrices and the one-loop masses. Effective pro- 
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duces the tree-level matrices and the corresponding diagonal form in order to perform 
all its other calculations. The full one-loop correction requires performing an additional 
diagonalization for each set of parameters one would want to compute the masses at. 
Since the diagonalization of the matrices is one of the most CPU intensive calculations 
performed by Effective, doing the full one-loop corrections is very time consuming. 

6.5.3 Effective Potential 

At this point the model has been completely defined and the tree-level masses have been 
computed. It is now possible to find the effective potential and minimize it over a set 
of parameters. The analytic tree-level effective potential is generated by Effective and can 
be returned. Effective also provides a function ex extremizePotential (vector<numeric*>) 
that goes through the parameters given in the input vector and minimizes the potential 
over these parameters using the Powell [11] minimization routine. The resulting mini- 
mum potential is returned and the values of the input parameters are set to the values 
that give a minimum. For example the tree-level effective potential of the electroweak 
model is given in figure 11.31 Passing the VEV parameter v into this routine returns 
— 1.69 X 10^ GeV^ for the potential and sets the value of v to 246.0 GeV. 

The one-loop potential can also be minimized in Effective. This uses the last calcu- 
lation of the mass matrix (whether it is tree-level or one-loop) to compute the one-loop 
correction given in ()6.38|1 . This can then be minimized in the same manner as the 
tree-level potential. 

Figure iniZI shows a plot of the MSSM potential from Effective with the given parameter 
set {/i = 30.0 GeV,m|^^ = 100.0 GeV^m^^ = -609.0 GeV^ 6 = 7000.0 GeV^} and the 
two gauge group couplings are set from the SM. In this plot we set vf + V2 = = 
(246.0 GeV)^ and tan/5 = ^. The tree-level effective potential of this model is derived 
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Figure 6.7: The effective potential of the MSSM for the parameter set {// = 
30.0 GeV,m|^^ = 100.0 GeV^m|^^ = -609.0 GeV^6 = 7000.0 GeV^} plotted in terms 
of (5. This potential can be minimized in Effective over the parameters, vi and V2. 

in Effective and is given by 

K. = ^ (/' + 'A) („t + „S - ^ + ^ + + !M + !|!. (6.90) 

From this equation we can see with the right set of parameters, the h term will push the 
potential below zero and require a particular combination of the VEVs to minimize. 

6.5.4 Renormalization Group Equations 

The last feature of a model that is generated automatically by Effective is the RGEs. 
Each parameter of the model can have its default value set at its own scale. The RGEs 
then evolve these parameters so that the values are all used consistently at any desired 
scale. 

As was shown earlier the (5 functions for each parameter are generated using the 
one-loop effective potential and the two-, three- and four-point Green's functions and 
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the anomalous dimension of the external fields. 

The solution to the boundary value problem can be given as follows. Start with a 
scale X — Mj. The boundary conditions at this scale can be given by 

f{p{x)) = 0. (6.91) 

where p is the set of parameters which are defined at scale x. The boundary condition 
at a new scale, z, is 

g{q{z)) = 0, (6.92) 

and again, q is the set of parameters defined at scale z. The full set of parameters of the 
model are given by {p{t), q{t)} where p{t) is constrained at t — x and q{t) is constrained 
at t — z. The parameters at the two scales are related by two functions, IF and i?^ by 

q{z) = R'^{q{x),p{x)) 

p{z) = RPiq{x),p{x)), (6.93) 

where these functions are determined by the RGEs. Using these equations we can 
numerically solve for the different parameters. 

There are three numerical methods that have been considered in Effective [114]. The 
first is known as the "shooting" method [11]. Using this method one computes g{q{z)) 
numerically by use of the function R*^. This amounts to numerically finding the solution 
to 

g{R\q{x),p{x)))^Q, (6.94) 

where p{x) is the boundary values of the parameters at scale x. This just requires solving 
this equation over the set of parameters q{x). 

The "shooting" method cannot be used in SUSY, however, as SUSY requires the 
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input to be given at three scales: Mz,Ms and M^. Another reason that this method 
cannot be used for SUSY theories is that the boundary conditions of a SUSY theory 
cannot be conveniently expressed as is done in ()fj.92|l . 

An alternative method [11] is the "drift" method. This method is a modified version 
of the "shooting" method. Here a guess qx = q{x) is taken and p{x) is determined from 
()(i.91|l . q{z) and p{z) are then given from ()fi.93j) . The boundary conditions at scale z 
are imposed to give q{z). These are then run back up to the scale x by the inverse of 
fl(i.93|) . This gives a new set of values, q'{x) and p'{x). The conditions at scale x are 
then imposed again. This process defines a recurrence relation q'{x) = Q{qx) of which 
the desired values of p and q will be a fixed point. The problem with this method is that 
the physical values might not be the only fixed point and there is no guarantee that this 
fixed point is stable. 

A third method, which has been implemented in Effective, handles the case where 
the fixed point is not stable. In this case we can instead numerically solve the equation 

Qiq.) -q. = 0, (6.95) 

using the Newton- Raphson [11] method. This method can be rather slow but if one 
makes reasonable guesses for some of the less important parameters the parameter space 
can be reduced and convergence is quicker. 

We have already seen the running of the gauge group couplings in fig. l6.5l produced by 
Effective. These were generated by setting the fields and groups of the model. Effective 
then computes the contribution to the (3 function from each field of the model. The 
coupling is evolved according to the (3 function generated. 

Figure 16.81 shows the evolution of the b mass given by the f3 function for the Yukawa 
coupling in the SM. This has been evolved using /3 functions produced by Effective and 
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Figure 6.8: The running of the b mass. This is given by the (3 function of the Yukawa 
couphng. 

the Runge-Kutta [11] method of differential equation evolution. The initial condition 
for this figure is mb{mz) = 3.0 GeV. Again, the /3 function is generated by Effective. 
Each of the anomalous dimensions of the Yukawa coupling are added together. This is 
the sum from ()6.81|) . The B factors have been calculated from the IPI contributions of 
the vertices given in fig. 16.61 

6.5.5 Future Extensions of Effective 

Before we discuss the extensions of Effective we first summarize what it can already do 
and how a user could use Effective without any extensions. 

Effective is able to generate the full Lagrangian of a SUSY model based on the groups, 
fields, superpotential and SUSY breaking terms. This Lagrangian is used to generate 
the effective potential, mass matrices and RGEs of the model. At this stage a user can 
input their parameters over a range of scales and evolve them to the same scale. The 
mass spectrum at tree- level or one- loop can be evaluated given the parameter set. The 
user could then automate the study of the mass spectrum over a range of parameters. 
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The code of Effective allows for several extensions. The most useful one in producing 
physical results is allowing additional diagrams to be coded. The library provides a class 
Diagram which provides several useful functions. This class allows a user to pass in the 
external fields of a diagram and a set for each of the internal propagators. The class 
provides functions which find the couplings of a set of fields given as either interaction 
states, mass states or combinations of the two. This is done by computing the coupling 
of mass states based on their mixing matrices. The couplings can then be used in the 
computation of a diagram by inheriting the Diagram class and implementing the virtual 
function ex functionO function. The diagram is then evaluated by iterating over all 
the internal propagator fields and summing the result of each diagram by a call to the 
ex evaluate function. Implementing ones own diagrams would allow Effective to 
be used to produce results of physical calculations (like cross-sections) over a range of 
parameters. 

The implementation of a users own diagrams isn't the only possible extension. The 
classes GaugeField and MatterField implement the set of expressions given in ()6.861 
I6.87|) . If someone wanted to consider a model which needed other terms (such as 
N ^ 1 SUSY) they would need to inherit these classes and override the function ex 
interactionO . They could then add new terms that depend on the fields. Unfortu- 
nately, this isn't all that would be needed for 7^ 1 SUSY. There would be a need to 
inherit the Model class and change several functions. But this is still a feasible task. 

One may also consider extra dimensional models. These too could be implemented in 
Effective. One would have to implement a mechanism for producing a four dimensional 
effective Lagrangian from the input (groups, fields, superpotential, plus any extra inter- 
actions) but this too would be a feasible task. It may be that an implementation such 
as this could only be done for a specific type of compactification scheme. But would still 
be useful for studying both SUSY and non-SUSY models derived from an effective four 
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dimensional Lagrangian. 

Another extremely useful extension that would be desirable is implementing a method 
of deriving the Feynman rules from the Lagrangian. All of the elements are in place to 
do such a task but would require quite a lot of work. These rules could be output in 
a format suitable for reading into an matrix element generator, like FormCalc, or one 
could even implement a matrix element generator as another extension of Effective. This 
would be a very complicated project, however. 

Effective in its current form can mainly be used to study the mass spectrum of a 
model once the VEVs that minimize the effective potential are determined. This mass 
spectrum depends on the input parameter set which can be given at several scales and 
the RGEs of the model can be used to match them at the scale the mass spectrum is 
desired at. Future extensions of the software would allow this mass spectrum to be used 
in calculations of cross-sections and other properties. More substantial extensions could 
allow the software to be used on 7^ 1 SUSY and extra dimensional models. Finally the 
development of Feynman rules for a model could be implemented with an eye towards 
use in matrix element generators. 



Conclusion 



This thesis has presented two new software packages, Herwig++ and Effective. Herwig++ 
is a Monte Carlo event generator. In its current state, this generator is able to generate 
e+e^-annihilation events. Chapter El has developed a new set of variables for the parton 
shower and chapter IHl has presented an improved hadronization model. A full discus- 
sion of the implementation specific features of Herwig++ have been given in chapter 0] 
Results of the new models and the package as a whole have been shown in chapter El 
These results give us great confidence that the new shower and hadronization are a 
better description of the physics. 

Currently, we are working on Herwig++ to generate initial-state showers. This will 
enable us to describe hadron-hadron and lepton-hadron events. The aim is to have this 
fully functional for use with Tevatron and LHC studies. Once the initial-state shower is 
fully functional a new underlying event model can be implemented as well as improving 
the hadronic decays. Work is already underway to include spin correlations into the 
shower and decays. New matrix elements for various processes can be implemented and 
SUSY particles and processes can be included. As the history of HERWIG has shown, 
there will continue to be a lot of potential for future research in Herwig++. 

Effective is a model building program. This program has been designed to work with 
= 1 SUSY models but future extensions could move it beyond just these models. This 
program automates the process of determining the mass spectrum of a model for a given 
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parameter set. The mass spectrum generated with Effective can be read into Herwig++ 
to define the masses of the supersymmetric particles in Herwig++. This would provide 
an automated way of studying SUSY at colliders over a range of parameters. Currently, 
only special purpose programs, hke SOFTSUSY [116], are able to do something similar. 
Future extensions of Effective may provide more functionality and make it an extremely 
powerful and useful tool for studying models. 



Appendix A 



Herwig++ 



This appendix contains a few relevant issues pertaining to the use and future develop- 
ment of Herwig-|— 1-. I first give an example analysis program to count the number of 
7r° in an event. Following that is a brief discussion of how to change the parameters in 
Herwig-|-+. Finally, 1 give the skeleton structure of how one would implement one's own 
matrix element. 



A.l Counting Pions 

This section provides a description of the functions needed to run a simple Herwig-I— |- 
program. The program given here will simply count the number of neutral pions in the 
final state. 

Herwig-|— I- can be run inside an exception handling clause. In C-|— I- this is the state- 
ment 

try { ... your code here . . . } 

catch(...) { // do something with the exception } 
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Herwig++ and ThePEG will throw exceptions to explain why the program has termi- 
nated. If these aren't caught the message associated with the end of the program won't 
be known and the reason for the termination will remain a mystery. Therefore, the main 
program to run Herwig-|— |- will be enclosed in such a statement and the message from 
the exception will be printed. The exact form we use is 
try { ...generate events /analyze ...} 
catchCstd: : exception & e) { 
cerr « e.whatO « endl; 
return 1; 
} catchC . .) { 

cerr « "Unknown exception\n" ; 
return 2; 

} 

We now move on to discussing what is put in the "generate events/analyze" sec- 
tion. A helper class HerwigRun has been developed for Herwig-I— 1-. This class takes the 
command line arguments (int argc, char **argv) and sets up Herwig++ to either 
initialize, read or run. These different stages are discussed in the next section. The 
command to create a HerwigRun class is 

Herwig: : HerwigRun hw(argc,argv) ; 
Once this class is created generating an event is straightforward. We first must check 
the status of the program, given by the functions isRunModeO and preparedToRunO. 
Wc can then iterate over the number of events to generate given by the function getN() . 
To generate an event we simply call generateEvent (). The code to do all of this is 
if (hw . isRunMode ( ) && hw . preparedToRun ( ) ) { 
forCint i = 0; i<hw.getN(); i++) { 
hw . generateEvent ( ) ; 
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At this point we have now generated an event. This is where the analysis code must 
be implemented. The HerwigRun class also offers a function to retrieve the particles from 
the last event generated. ThePEG: :tPVector getFinalStateCint i) is this function. 
If i is not given (e.g. getFinalStateO) then the final state particles of the event are 
returned. If the i argument is given this function returns the particles in the final state 
given by the step i. 

Returning to our example of counting the number of ir^ particles in the final state 
we pass the vector of final state particles to a counting function which we define, int 
countPions (ThePEG: :tPVector particles). This function is fairly straightforward 
and is given by 

int countPions (ThePEG: :tPVector particles) { 
ThePEG: :tPVector: : iterator it; 
int count = 0; 

f or(it=particles .beginO ; it ! =particles . end() ; it++) { 
if ((*it)->id() == ThePEG: :ParticleID: :piO) count++; 

} 

return count; 

} 

This function has rehed on the class ThePEG: : Particle which defines several func- 
tions of the particle. Here we have used the long idO function to determine the PDG 
code of the particle. There are many more functions that would be useful for analysis. 
The user is encouraged to read the documentation of the Pcirticle class for more details. 

The HerwigRun class also defines functions which return some of the elements from 
ThePEG. These are things like the Step or the CollisionHandler from which the user 
is able to obtain all the information of the event. 
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The repository is a feature of ThePEG that allows the parameters of the program to 
be changed without having to recompile the code. These parameter files can then be 
exchanged within collaborations to ensure that researchers are working with the same 
parameters. 

The repository can be thought of as the input files to Herwig++. These are writ- 
ten in a human readable form, though the repository has a simple syntax of its own. 
The input files can then be saved in machine readable form for faster access on sub- 
sequent runs. In Herwig-|— |- there is a default input file, HerwigDef aults . in. This 
sets up all the relevant objects and instructs the program to allocate memory for 
the various parts of the event generator. This initial setup can be saved in a ma- 
chine readable file, HerwigDef aults .rpo. The HerwigRun class reads in the command 
line argument init. This instructs the class to create the .rpo from the .in file. 
HerwigDef aults . in also reads Shower, in, Hadronization. in and Decays, in as well 
as the particles in ThePEGParticles . in and HerwigParticles . in and the decay modes 
from HwDecays . in. 

The HerwigDef aults . in currently defines two types of generators. One for LHC- 
like events and another for LEP-like events. These two generators can then have their 
parameters changed (such as cm. energy or various cuts) by editing the files LHC . in 
or LEP.in, respectively. These modifications arc read in by the read command line 
argument. This reads in the appropriate .in file and produces a .run file. The .run is 
checked to ensure that all the relevant objects and links have been set so that the event 
generator can run. 

The last stage is to actually run some events. Using the command line argument 
run reads in the appropriate .run file and generates events. There is a parameter 
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NumberOf Events in each generator object that sets the maximum number of events to 
run. The HerwigRun class provides for a new number of events to be given at runtime. 
This number must be less than or equal to the parameter given to the NumberOf Events 
field. To instruct the program to use a new number simply use the command line 
argument -N # when starting the program. 

To summarize there are three ways in which a program that uses a HerwigRun class 
can by used. These are given as 
Prograin_name init 

Prograiii_name read Generator_f ile . in 

Prograin_name run Generator_f ile . run -N # 
There are more commands that can be given in the command line but these mostly 
change some of the lower level instructions. The full set of commands is 

Prograin_name init I read I run [-N num-events] [-seed random-generator-seed] 
[-d debug-level] [-dHw herwig-debug-level] [-1 load-path] [-L first-load-path] 
[-r repo-file] [-i initialization file] [run-file] 

To change the value of a parameter you find the line in the input file that this is 
governed by. For example if we wanted to turn the initial-state radiation off we would 
look in the Shower . in file. This file has the line 

set theSplittingGenerator:OnOff ISRMode 1 
If we simply edit this file and change the 1 into a this will turn the initial-state radiation 
off. Most of the commands are straightforward like this. There are other commands 
in the repository than set. A relevant subset is create, set, mkdir,cd, library. The 
create command is used to define a new object given by a class. The set command is 
used to change one of the parameters of that object. The repository is structured like a 
file system. You can create the objects you want in directories in order to keep things 
ordered. The mkdir command creates a new directory and the cd command changes 
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to the given directory. The last command is the library command. This is used to 
dynamically load a library. This must be done to define objects that are in a library 
that hasn't been loaded yet. For example, if we want to define an Amegiclnterf ace 
object we would have to load the libHwAmegic. so library first. 



A. 3 Matrix Element Development 

In this last section I show how a user would implement their own matrix element. The 
class MEBase defines the lowest level of matrix element abstraction. To implement a 
matrix element one must at least inherit this class and define the functions 

unsigned int order InAlphaSO const 

unsigned int orderlnAlphaEWO const 

double me2() const 

Energy2 scale () const 

generateKinematics (const double r) 

CrossSection dSigHatDRO const 

void getDiagramsO const 

Selector<const ColourLines*> colourGeometriesCtcDiagPtr diag) const 

This list looks quite daunting. ThePEG has already generated these functions for special 
types of matrix elements. For example the class ME2to2Base as able to define the 
generateKinematics (), dSigHatDRO as well as the scale () for all 2 ^ 2 processes. A 
further subclass ME2to2QCD defines order InAlphaSO and orderlnAlphaEWO functions 
as well as providing routines for general functions like determining the number of active 
fiavours at a given scale. 

If we take an example to inherit the ME2to2QCD class this means we only need to 
define the functions 
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double me2() const 

void getDiagramsO const 

Selector<const ColourLines*> colourGeometriesCtcDiagPtr diag) const 

The function double me2() const is self explanatory. This is just the value of the 
matrix element for the set of currently generated points. To access the set of points there 
is a function meMomentaO which returns a vector of LorentzSMomentum objects. This 
vector contains the momentum generated for all the incoming and outgoing particles. 

The getDiagramsO const class is used to indicate to the matrix element what 
are the diagrams used to produce the matrix element (squared). There is a function 
add(DiagPtr) which adds a diagram to the list. There is a helper method defined to 
define a diagram. For example 

new_ptr ( (Tree2toNDiagrain(2) , q, qb, 1, gamma, 3, 1, 3, lb, -1)) 
define a 2 ^ 2 process where the q and the qb are the incoming particles. The 1 indicates 
that the gamma particle is the child of the first incoming particle. The 3 , 1 , 3 , lb series 
indicates that both the 1 and lb particles are the children of the gamma particle. The 
last argument of -1 indicates that this is the end of the diagram definition. This is 
a fiag that must be negative but can have different values. The different values are 
used to determine the different types of diagrams. If one wanted to weight the different 
diagrams they would need to overload the Selector< Diagramlndex> diagrams (const 
DiagramVector &) const function. In this function one can iterate over the diagrams 
and add them into the Selector object with a weight. When a diagram is selected the 
weights are taken into account. 

The last function that must be implemented is the Selector<const ColourLines 
*> ColourGeometriesCtcDiagPtr diag) const function. This function allows you to 

set the different colour connections. These can also be given a weight which the particular 
colour connection implemented is chosen from. 
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Having implemented all of these functions the matrix element is defined and can be 
used in Herwig++. The typical place to use the matrix element is in the hard subprocess. 
It is possible to use a matrix element in other processes, such as in the decays, but that 
requires more development and working with the code on a lower level. 



Appendix B 



Effective 

This appendix is devoted to coding issues of Effective. I give here an example of coding 
the MSSM model in Effective and how one can use the model to get the mass of the 
neutrahnos. I also show how to evolve the RGEs to provide the scale dependence of the 
parameters. At the end I give a skeleton outline of how one would implement one's own 
groups and diagrams. 

B.l MSSM 

The main driver of Effective is the class Model. This class is what directs all the calcula- 
tions and provides all the components. This class that must be inherited to implement 
ones own model. The methods that need to be defined are 

void createGaugeGroupsO 

void createGaugeFieldsO 

void createMatterFieldsO 

void addOtherTermsO 

ex superPotentialO 
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The createGaugeGroupsO function specifies what groups are used in the model. 

For example in MSSM and the SM we have only the SU{3)c, SU{2)l and U{1)y groups. 
These can be implemented via the function 

void MSSM: : createGaugeGroupsO { 

addGaugeGroupCnew UlGroupC'Ul" , "{g'}", this, UlY)); 
addGaugeGroupCnew SU2Group("SU2" , "{g_W}", this, SU2w)); 
addGaugeGroupCnew SU3Group("SU3" , "{g_3}", this, SUSc)); 

} 

The objects UlY, SU2w and SUSc are all integers that define the 'line' of the group. The 
two strings for each group define the name of the group coupling in the standard text 
output and the I^TgX output, respectively. The groups can then be accessed at a later 
point using the standard text string (e.g. getGroupC'Ul")). 

The next function that needs to be implemented is createGaugeFieldsO. This 
function defines the fields that mediate the interactions. In the MSSM these are the 
superpartners {B'^,B), (Wf , Wj) and {A'^,Aa). These are given by 

void MSSM: : createGaugeFieldsO { 

addPieldC'B" ,new GaugeFieldC'B" , "B" ,Utils: :LorentzVector, 

getGaugeGroup ( "Ul " ) ) ) ; 
addFieldC'Bino" ,new GaugeFieldC'Bino" , "\\tilde{B}" ,Utils : :WeylSpinor , 

getGaugeGroup ("Ul") ,getGaugeField("B") ) ) ; 
addField("W" ,new GaugeField("W" , "W" ,Utils: :LorentzVector, 

getGaugeGroup ( " SU2 " ) ) ) ; 
addField("Wino" ,new GaugeField("Wino" , "\\tilde{W}" ,Utils: :WeylSpinor, 

getGaugeGroup ("SU2") ,getGaugeField("W") ) ) ; 
addField("gluon" ,new GaugeField("A" , "A" ,Utils: :LorentzVector, 
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getGaugeGroup ( "SU3" ) ) ) ; 
addFieldC'gluino" ,new GaugeFieldC'Aino" , "\\tilde{A}" ,Utils : : WeylSpinor, 

getGaugeGroup ("SU3") ,getGaugeField("gluon") ) ) ; 

} 

The addPieldO function takes a string as an index for the field and a pointer to 
a field. This same function is used to add the MatterFields. Creating a GaugeField 
requires two strings. Again the first is the text output string and the second is the 
output. These strings are followed by the spin of the particle. This is actually 
the (2s + 1)(— l)^'* factor. This has been conveniently defined in the class Utils for the 
scalar, fermion and vector spin types. The spin is then followed by a pointer to the 
gauge group that the field is the mediator of. The last argument, which is optional, is 
the SUSY partner of the field. We can see above that the SUSY partner is set for the 
fermion but not the vector. This is because one field must be created first. Once created 
the partner can be set. Setting it in the creation of the fermion also sets it in the vector, 
so this syntax sets both partners to each other with only one reference to a partner. 

After the addGaugeFields function is defined we implement the addMatterFields 
function. This function sets all of the remaining fields in the model. Again looking at 
our example from MSSM we can write the function as 

void MSSM: :createMatterFields() { 
numeric half (1,2); 
numeric sixth(l,6); 
numeric third (1,3); 
numeric twothirds(2,3) ; 
// Adding leptons/sleptons 

addFieldC'leptonL" ,new MatterFieldC'L" , "Well" , Utils: : WeylSpinor, famsize, 

getGaugeGroup ("Ul") ,-half , 
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getGaugeGroup ( " SU2 " ) , 1 ) ) ; 
addFieldC'sleptonL" ,new MatterFieldC'sL" , "\\tilde-[\\ell}" , 

Utils: : Scalar, famsize, 
getGaugeGroup ( "Ul " ) , -half , 
getGaugeGroup ( " SU2 " ) , 1 , 
getMatterFieldC'leptonL"))) ; 

addPieldC'leptonR" ,new MatterFieldC'eR" , "e_R" , Utils: rWeylSpinor, famsize, 

getGaugeGroup ("Ul") ,-1)) ; 

addField("sleptonR" ,new MatterField("seR" , "\\tilde{e_R}" , 

Utils : : Scalar , famsize , 
getGaugeGroup("Ul") ,-1, 
getMatterField("leptonR"))) ; 

// Adding quarks and squarks 

addField("quarkL" ,new MatterField("Q" , "Q" , Utils: :WeylSpinor, famsize, 

getGaugeGroup ("Ul") .sixth, 
getGaugeGroup ( " SU2 " ) , 1 , 
getGaugeGroup ("SU3") , 1)) ; 
addField("squarkL",new MatterField("sQ" , "\\tilde-[Q}" , 

Utils: : Scalar, famsize, 
getGaugeGroup("Ul") , sixth, 
getGaugeGroup ( " SU2 " ) , 1 , 
getGaugeGroup ( " SU3 " ) , 1 , 
getMatterField("quarkL"))) ; 



addField("uR" ,new MatterField("uR" , "u_R" , Utils: :WeylSpinor, famsize. 
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getGaugeGroup ( "Ul " ) , twothirds , 
getGaugeGroup("SU3") ,1)) ; 
addFieldC'suR" ,new MatterFieldC'suR" , "\\tilde{u_R}" ,Utils: : Scalar, famsize, 

getGaugeGroup ( "Ul " ) , twothirds , 

getGaugeGroup ( " SU3 " ) , 1 , 

getMatterFieldC'uR"))) ; 

addPieldC'dR" ,new MatterFieldC'dR" , "d_R" ,Utils: :WeylSpinor, famsize, 

getGaugeGroup ("Ul") , -third, 

getGaugeGroup ( " SU3 " ) , 1 ) ) ; 
addField ( " sdR" , new MatterField ( " sdR" , " \\t ilde{d_R} " , Ut ils : : Scalar , famsize , 

getGaugeGroup ("Ul") , -third, 

getGaugeGroup ( " SU3 " ) , 1 , 

getMatterField("dR"))) ; 

// Higgs fields 1 is upper, 2 is lower 

addField ("HI", new MatterField("Hl" , "{H"l}" ,Utils : :Scalar,l, 

getGaugeGroup ( "Ul " ) , -half , 
getGaugeGroup ("SU2") ,1)) ; 
addField ("sHl", new MatterField("sHl" , "\\tilde{H"l}" ,Utils : rWeylSpinor , 1 , 

getGaugeGroup ("Ul") ,-half , 
getGaugeGroup ( " SU2 " ) , 1 , 
getMatterField("Hl"))) ; 



addField ("H2", new MatterField("H2" , "{H~2}" ,Utils : : Scalar,!, 

getGaugeGroup ( "Ul " ) , half , 
getGaugeGroup ("SU2") ,1)) ; 
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addField("sH2", new MatterField("H2" , "\\tilde{H"2}" ,Utils : : WeylSpinor , 1 , 

getGaugeGroupC'Ul") ,half , 
getGaugeGroup ( " SU2 " ) , 1 , 
getMatterField ( "H2 " ) ) ) ; 

// Now add higgs Vev 

Parameter upsilonl = addParameterC'HiggsVevl" , "\\upsilon_l" ,220.0, 

Parameter: :vev) ; 

Parameter upsilon2 = addParameter("HiggsVev2" , "\\upsilon_2" ,220.0, 

Parameter: :vev) ; 

addVev ( "HiggsVevl " , getPield ( "HI " ) , 1st (get Index ( "HI " , " SU2 " ) ==1) , 

(upsilonl+Model : : star) ) ; 
addVev("HiggsVev2" ,getField("H2") , 1st (getIndex("H2" , "SU2")==2) , 

(upsilon2+Model : : star) ) ; 
addVevParameter (upsilonl) ; 
addVevParameter(upsilon2) ; 

The meaning of this code is easy to see. Again we see the last (optional) argument 
in creating a MatterField is the superpartner. We can see now that there are some 
new arguments to these fields though. There is a f amsize argument. This is an integer 
which defines how many generations the field has. For example in MSSM we would set 
f amsize to 3. Lastly we see that after each gauge group is a number. This is the 'charge' 
of the group. Most groups just have a charge of 1. The U{1)y does not, however. This 
has fractional charges. The class numeric is one from GiNaC which defines rational 
numbers, without rounding errors. These expressions will be printed in as they 

are (e.g. half = |). 

In the last part of the code we see that we have created the VEVs. The class 
Parameter is used in Effective to represent all of the parameters of the model. The 
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Parameter class is printed out in symbolic form during printing but is evaluated to its 
value during calls to evalf (). The object Model: :star is a placeholder which is where 
the relevant field is placed. If we wanted to also add an imaginary VEV would could 
use a similar syntax imagVEV+Model : : star and in this instance Model : : star is the 
imaginary part of the field. 

As discussed in the text the parameters could be added in a different way. If we 
wanted to work with v and tan /? rather than vi and V2 we could have used the code 

Parameter upsilon = addParameterC'HiggsVev" , "$\backslash\backslash$upsilon" ,246.0, 

Parameter: :vev) ; 

Parameter beta = addParameter ("beta" , "$\backslash\backslash$beta" , 1 . 1 , 

Parameter: :vev) ; 

addVev ( "HiggsVevl " , getField ( "HI " ) , 1st (get Index ( "HI " , " SU2 " ) ==1) , 

(upsilon*sin(beta)+Model : :star)) ; 
addVev("HiggsVev2" ,getField("H2") , lst(getIndex("H2" , "SU2")==2) , 

(upsilon*cos(beta)+Model : :star)) ; 
addVevParameter (upsilon) ; 
addVevParameter (beta) ; 

but as discussed in the text, this type of setup will not properly define the RGEs and 
the minimization routine may not succeed. 

In the MSSM there are many terms which must be added from both the superpo- 
tential and from the soft breaking terms. We won't give them all here but instead will 
present an example of how to add a term. We will look at how to add the slepton mass 
term. This is given by 

ex mOLL = addFamilyMatrix("mOLL" , "{m'^2_L}", f amsize) ; 
vector<idx*> sumlndex; 
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idx i = Utils: :f amilylndexCOjfamsize) ; 

idx j = Utils: :f amilylndexCl ,f amsize) ; 

MatterField *sleptonL = getMatterFieldC'sleptonL") ; 

ex li = sleptonL->expression() ; 

ex Ij = Utils :: conjugate (li) ; 

Ij = Ij . subs(sleptonL->f a]nilylndex()==j) ; 

sumlndex = sleptonL->getIndices() ; 

sumlndex . push_back (& j ) ; 

ex sleptonLTerm = Utils: : real (Utils : :suinIndices(mOLL*li*lj , sumlndex) ) ; 
add (-slept onLTerm) ; 

At first this looks quite horrific. But oucc the different parts are explained it is 
obvious what everything means. The first term is used to define the slepton mass term. 
Since there are three slepton generations the mass-squared coupling must be defined as 
a 3 X 3 matrix of parameters. The addFamilyMatrixO function provides this. This 
function creates a matrix whose elements are parameters. The parameters are created 
with the indices appended to the name (e.g. the mj^n element is given by the parameter 
mOLLll). When the matrix is created it is created as an identity matrix. 

The next line defines a vector in which all the indices which are summed over are 
placed. This is followed by defining two family indices. By default the FamilyMatrix 
class uses the first two indices from the Utils: :f amilylndexO function. The Utils 
class defines some universal symbols for different indices. A call to the f amilylndex 
function returns one of the symbols (given by the first integer) as an index where the 
dimension of the index is given by the second argument (e.g. the number of families). 
As a note, there is a integer Utils: :max_indices which defines how many symbolic 
indices are available. 

The definition of the two family indices, i and j, is followed by a retrieval of the 
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MatterField pointer to the slepton. If we refer back to our definition of the slepton field 
we can see that the first string is "slepton" . This is the tag that is used to retrieve the field 
at a later time. There is also a function getGaugeField which behaves in a similar man- 
ner, only it returns a GaugeField pointer. After the retrieval of the slepton pointer we 
now get the expression of the slepton. As Effective seperates the real and imaginary parts 
of the scalars and fermions, the expression is ■^{TZ{f)+iT{f)). Also the function expres- 
sion returns the fields with each index attached. For example a left-handed quark would 
have its family index as well as the SU (2) ^ and SU{3)c indices added. A vector also has 
a Lorentzlndex attached. The function Utils: : conjugate then replaces all i's with 
—i. The default family index of a field is taken from Utils : :f ainilylndex(0,f amsize). 
Since we want to contract the family indices of the mass-squared coupling we must re- 
place this index with Utils: : familylndexd ,f amsize) in the conjugate expression. 
This is done via the subsO command. This command takes two forms. The first is to 
use the == sign. This sets the argument on the left to be replaced by the argument on 
the right. This is useful when making one substitution. When many substitutions are 
desired we can instead use the form subs (1st ,1st). The first 1st is replaced by the 
second. A 1st is a class from GiNaC and is just a sequence of expressions. 

After the expression for and 0* have been retrieved we want to contract all the 
indices and create the expression m|,j^-0*0i. To do this we store all of the indices given 
in the slepton into the sumlndices vector we created earlier. We also realize that this 
vector doesn't contain the extra family index in m|, and now in Ij, so we also add the 
index to the back of the list. Once the sumlndices vector contains all the indices which 
we are summing over we can call the routine Utils: : sumlndices (ex, vector<idx*>). 
This routine iterates the list of indices over all permutations and evaluates the ex at 
each set of index values. The sum of all these permutations and evaluations is returned. 
The last command on this line is the Utils: :real. This returns all the parts of the 
expression that are not multiplied by i. In this example we shouldn't need to call this 
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command as 0*0 should be a real quantity. It is also worth noting that instead of coding 
an expression as exp + c.c. it is easier and faster to write 2*Utils: :real(). 

In the end we call the addO function. This adds the expression to the Lagrangian. 
If we were developing terms of the superpotential we would instead sum all the contri- 
butions from all the terms and return it. We don't want to call the add() routine with 
the superpotential. 

B.2 Plotting Effective Potential and Running Cou- 
plings of MSSM 

Once we have developed the class MSSM with all the relevant groups, fields and terms, 
we can work with this program. The following code is used to create the model and 
initialize it 

MSSM mssm; 

Model: :readCommandLine(argc,argv,&mssm) ; 
mssm. couplings ("mssm. gar") ; 
mssm. initialize ; 

The first line creates an MSSM object. The second line calls a static function from the 
Model class which reads in standard command line arguments and sets the relevant 
parameters of the model. The possible arguments are 

-p # - Sets the evaluation of the effective potential to tree-level (0) 
or one-loop (1) 

-t # - Sets the evaluation of the tadpoles to tree-level (0) or 
one-loop (1) 
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-a # - sets the evaluation of the mass matrices to tree-level (0) or 

one-loop (1) 
-r - Rederives the couplings from scratch 
-s - Saves the couplings to the file provided 
-h - Prints a help message 

The third line of code provides the model with a file to read the couplings in from 
and save to, depending on what the program has been specified to do. The last line 
initializes the model. This means generate (or read) the couplings, find all the tree level 
mass matrices and generate the RGEs. At this stage the model is able to be used. 

At this point we will want to read our parameters in from a file. This is done by 
the loadParametersO function. The code that would read the parameters from the file 
MSSM.dat is 

ifstream parin; 
par in. open ("MSSM.dat") ; 
mssm.loadParameters (parin) ; 
parin. closeO ; 

This file has a simple structure. It is tab delimited where the first entry is the name 
of the parameter (the text name given during its definition) and the value of the pa- 
rameter. There is also a parameter defined for every model. This is renormScale and 
defines the current renormalization scale. It is defaulted to 91.2 GeV (the Z° mass) and 
when a parameter is read in it is assumed to be entered at the scale currently given by 
renormScale. So if one wanted to input parameter a at 91.2 GeV but parameter b at 
200 GeV they would use 



a 50.0 
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renormScale 200.0 
b 20.0 
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Once the parameters have been read in one may want to evolve the RGEs to a scale. 
At first the boundary conditions must be determined at one scale. This is done by call- 
ing mssm.rgesO .matchScalesO. This routine goes through the boundary conditions 
entered and uses the Newton-Raphson method to find a suitable set of parameters that 
match all the boundary conditions given. Once that is done the parameters can be 
evolved to the scale desired by calling mssm.rgesO .evolve (scale). 

Once the parameters have been determined at the scale desired the mass spectrum 
can be developed. This is done with a call to mssm . masses () . diagonalize () . This will 
diagonalize all the mass matrices at either tree-level or one-loop. There is an additional 
flag, Model : : CorrectMassStates which determines whether the one-loop corrections 
are applied directly to the mass eigenstates derived from the tree-level mass matrices, 
or if the one-loop corrections are applied based on the interaction eigenstates. If one 
wanted to get the mass of a particle they simply call mssm. masses () .massOf (field) 
where field is the expression with all the indices evaluated (e.g. instead of looking at the 
top quark mass, one would have to look at the left-handed quark, with SU (2) l index = 
1, family index = 3 and the SU{3)c index set to any number between 1 and 3). In the 
SM, for example, one would fully expect to find that the mass all of the colour indices 
of a particular quark are the same. But it would be possible to develop a theory where 
this is not the case. 

The last function that would be desired is to find the effective potential. This can 
be given as an analytic expression by a call to mssm.treePotentialO or at one-loop 
order by a call to mssm.potentialO (with the correct value of Model: :_potApprox 
set) . Using this we can use the code 
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of stream potPlot; 

potPlot. open ("MSSM_Pot.dat") ; 

double vl , v2; 

ex pot = mssm.treePotentialO ; 
double V = 246.0; 
double beta; 

forCbeta = -3.14; beta < 3.14; beta += 0.01) { 
mssm.getParameter ("HiggsVevl") = v*cos(beta); 
mssm.getParameter ("HiggsVev2") = v*sin(beta); 
potPlot « beta « "\t" « pot.evalfO « endl; 

} 

potPlot . closeO ; 
to generate fig. 16.71 



B.3 Entering a New Group 

When studying models different symmetries are desired. These are given by the group 
governing the interaction. Effective provides a way in which a new group can be defined. 
Inheriting the class GaugeGroup is how this is done. This class requires the user to 
implement 

constructor (string n, string cn, Model* m, char 1); 
ex structureConstant(Gaugelndex&, Gaugelndexfe, Gaugelndexfe) const; 
ex generator (Gaugelndexfe, Matter Indexfe, Matterlndexfe) const; 
Gaugelndex gaugelndex(int i=0) const; 
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Matterlndex matter Index (int i=0) const; 
bool isIndexO const; 
ex Cr(bool) const; 
ex C2(bool) const; 

The constructor defines the couphng (given by the text name n and the name 
cn). The char 1 is used to define the fine of the group. A constructor is generally of 
the form ClassName(n,cn,in,l) : GaugeGroup(n,cn,m,l){}. This just passes all the 
handling of these arguments to the GaugeGroup superclass. The functions C2(bool) and 
Cr(bool) are used to give the evalation of t^t^ — C2{r) ■ 1 and C(r) — Tr[t^t^]. If the 
boolean passed in is true then the adjoint representation is used. If it is false then the 
fundamental representation is used. 

Next one must define the indices. This usually amounts to defining a set of symbols 
for the adjoint and fundamental indices. The indices are then just idx classes with 
the symbol and the correct dimension set. The last remaining pieces are the structure 
constants and generators, given by structureConstant and generator. These require 
the definition of a tensor class from GiNaC. An example of an SU (3) structure constant 
is given here. The definition of the tensor is 

class SUSStructure : public tensor { 

GINAC_DECLARE_REGISTERED_CLASS (SUSStructure , tensor) 

public : 

void print (const print_context &c, unsigned level = 0) const; 
ex eval_indexed(const basic &i) const; 

}; 
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One can then define a helper function SU3_structure which returns the expression for 
the structure constant given the indices. This takes the form 



inline ex SU3_structure (const ex &il, const ex &;i2, const ex &i3) { 
if (!is_a<GaugeIndex>(il) I I ! is_a<GaugeIndex>(i2) || 
! is_a<GaugeIndex> (iS) ) 
throwCstd: : invalid_argument( 

"Indices to SU(3) structure must be of type Gauge Index") ) ; 

if (ex_to<GaugeIndex> (i2) . line () ! = 
ex_to<GaugeIndex>(i3) .lineO) 
throwCstd: : invalid_argument ( 

"Gaugelndices must be of the same line and type")); 

return indexed(SU3Structure() , sy_anti(), il, i2, 13); 

} 

The structure constant can then be defined by 

GINAC_IMPLEMENT_REGISTERED_CLASS (SU3Structure , tensor) 
DEFAULT.CTORS (SU3Structure) 
DEFAULT.ARCHIVING (SU3Structure) 
DEFAULT.COMPARE (SU3Structure) 

and the print and eval_indexed functions are 

void SU3Structure: : print (const print_context & c, unsigned level) const { 
if (is_a<print_tree>(c)) inherited: : print (c, level); 
else { 



APPENDIX B. EFFECTIVE 239 

ostreamfe os = c . s ; 
OS « "f"; 

} 

} 

ex SUSStructure :: eval_indexed (const basic &i) const { 
GINAC_ASSERT(is_a<indexed>(i)) ; 
GINAC_ASSERT(i.nops() == 4); 
GINAC_ASSERT(is_a<SU3Structure>(i.op(0))) ; 
GINAC_ASSERT(is_a<GaugeIndex>(i . op(l) ) ) ; 
GINAC_ASSERT(is_a<GaugeIndex>(i . op(2) ) ) ; 
GINAC.ASSERT (is_a<GaugeIndex> (i . op (3) ) ) ; 

const Gaugelndexfe gl = ex_to<GaugeIndex> (i . op(l) ) ; 
const Gaugelndex& g2 = ex_to<GaugeIndex>(i . op(2) ) ; 
const Gaugelndexfe gS = ex_to<GaugeIndex>(i . op(3) ) ; 

// Numeric Evaluation 

if (static_cast<const indexed &>(i) . all_index_values_are(inf o_f lags : : integer)) 
{ 

int i = ex_to<nuineric>(gl.get_value()) .to_int() ; 
int j = ex_to<nuineric>(g2.get_value()) .to_int() ; 
int k = ex_to<numeric>(g3.get_value()) .to_int() ; 
int ip = i; int jp = j; int kp = k; 
if(ip>jp) Utils: :swap(ip, jp) ; 
if(jp>kp) Utils :: swap (jp,kp) ; 
if(ip>jp) Utils: : swap (ip,jp) ; 
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static ex half = numeric (1, 2) ; 

double eye = Utils: :cyclicPermutation(i, j ,k) ; 

if(ip==l) { 

if(jp==2 && kp==3) return eye; 

else if(jp==4 && kp==7) return half*cyc; 

else if(jp==5 M kp==6) return -half*eye; 
} else if(ip==2) { 

if (jp==4 && kp==6) return half*eye; 

else if(jp==5 M kp==7) return half*eye; 
} else if(ip==3) { 

if (jp==4 && kp==5) return half*eye; 

else if(jp==6 M kp==7) return -half*eye; 
} else if(ip==4 && jp==5 && kp==8) return sqrt (numerie(3,4))*eye; 
else if(ip==6 && jp~7 && kp==8) return sqrt(numerie(3,4))*cye; 
return 0; 

} 

// No further simplification 
return i .holdO ; 

} 

This then completely defines the behaviour of the structure constants of the group. 
A similar thing can be done for the generators. 

B.4 Defining a Diagram 

Defining a diagram is a useful way to do computations with Effective. This is also 
straightforward. One simply inherits the Diagram class and implements the constructor 
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MyNewDiagramCDElementVec &v. Model *m) : Diagram (v,in) {} 
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and the functionO function. This function is used to evaluate the diagram at the 
current set of fields. The DElementVec is used to specify a set of fields for each external 
leg and propagator of the diagram. When the diagram is created the desired choice 
of fields for each part is given. For example we want to compute the cross section 
for e+e" —>■ qq at tree-level. We simply create some DiagramElements and put the 
appropriate fields in. 

DElementVec eeqq(4) ; 

exvector electrons, quarks; 

electrons . push_back (real_lef t_electron) ; 

electrons .push_back(imag_left_electron) ; 

electrons . push_back (real_r ight_electron) ; 

eeqqEO] = DiagramElement (electrons , true); 
eeqq[l] = DiagramElement (electrons , true); 
eeqq[2] = DiagramElement (quarks, true); 
eeqq[3] = DiagramElement (quarks, true); 

This creates the DElementVec object with all the electrons and all the quarks. This then 
gets given as input to the diagram ee2qq 

ee2qq myee2qq(eeqq,my_model) ; 

cout « "My matrix element is " « myee2qq. evaluate () « endl; 

The evaluate function goes over every permutation of the fields given in the four 
DiagramElements and computes a contribution to the matrix element. This contribu- 
tion is given by the function function. This must be defined so that each combination 
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of fields contributes the relevant part. This is done by determining the couphngs of the 
set of fields given. For example 

ex coupl23 = coupling(0, 1,2) ; 
ex coup234 = coupling (1,2, 3) ; 
ex coupl34 = coupling (0,2, 3) ; 
ex coupl234 = coupling (0, 1,2,3) ; 



The true statement in the definition of the DiagramElement indicates that the mass 
state couplings are to be used. If it were false the interaction couplings would be used. 
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